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Abstract. 

Topological states of fermionic matter can be induced by means of a suitably 
engineered dissipative dynamics. Dissipation then does not occur as a perturbation, 
but rather as the main resource for many-body dynamics, providing a targeted cooling 
into a topological phase starting from an arbitrary initial state. We explore the concept 
of topological order in this setting, developing and applying a general theoretical 
framework based on the system density matrix which replaces the wave function 
appropriate for the discussion of Hamiltonian ground-state physics. We identify key 
analogies and differences to the more conventional Hamiltonian scenario. Differences 
mainly arise from the fact that the properties of the spectrum and of the state of the 
system are not as tightly related as in a Hamiltonian context. We provide a symmetry- 
based topological classification of bulk steady states and identify the classes that are 
achievable by means of quasi-local dissipative processes driving into superfluid paired 
states. We also explore the fate of the bulk-edge correspondence in the dissipative 
setting, and demonstrate the emergence of Majorana edge modes. We illustrate our 
findings in one- and two-dimensional models that are experimentally realistic in the 
context of cold atoms. 
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1. Introduction 

Symmetries, their spontaneous breaking, and related order parameters were considered 
for a long time as the paradigm for understanding ordered states of matter. A paradigm 
shift was initiated in the late 1980s when the Landau-Ginzburg broken-symmetry theory 
of ordered phases — widely thought to be exhaustive — proved unable to characterize a 
new kind of phases with no local order parameter: topological phases, or phases with 
topological order [1, 2]. Instead of being distinguished by symmetries, topological phases 
are characterized by distinct values of a non-local, topological order parameter, and 
phases transitions occur whenever the topology changes, signaled by discontinuities 
in this topological invariant. The existence of topological order may be conditioned 
on the existence of symmetries. However, as long as topological order is present, the 
underlying system generally exhibits topological features, i.e., features that are robust 
against arbitrary (symmetry-preserving) quasi-local perturbations. 

Spectral gap and ground-state degeneracy are typical topological properties which 
have been theoretically shown to be robust for wide classes of Hamiltonians [3-6]. 
Whereas the spectral gap is a property of the bulk — as topological order itself — the 
ground-state degeneracy generally depends on the boundary conditions imposed at the 
edges of the system and on the existence of topological defects in the bulk (e.g., vortices 
in a supcrfluid). Most importantly, the degeneracy can be traced to the existence of 
zero-energy modes localized at the edges or bound to topological defects, which are 
robust topological features as well. These objects can exhibit exotic behavior under 
spatial exchange (or "braiding") such as non-Abelian statistics [7-9], which opens up 
exciting possibilities for practical applications such as topological quantum memories 
and topological quantum computation [9-11]. 

The search for topological phases exhibiting quasiparticlcs with non-Abelian 
statistics has brought p-wave paired superfluids and superconductors to the forefront 
of theoretical and experimental condensed-matter research [8,9,11-13]. Such systems 
have first been studied in two dimensions (2D), where they have been predicted to 
support topological phases with gapless chiral edges modes and quasiparticlcs known as 
Majorana zero modes, giving rise to Ising-type non-Abelian exchange statistics [8, 12, 14]. 
Following a seminal paper by Kitaev [15], the focus has moved more recently to networks 
of one-dimensional (ID) systems, which were shown to allow for similar topological 
features (non-Abelian statistics, in particular) as genuine 2D systems [16,17]. Recent 
proposals for solid-state [18, 19] and cold atom [20] systems have made it possible for 
Majorana zero modes to enter the experimental stage, with promising first results in 
solid-state devices [21-24] and the perspective of increased future experimental efforts. 

In recent years, the quest for topological states was extended to non-equilibrium 
systems, going beyond the Hamiltonian ground-state scenario. A first step in this 
direction was taken with periodically driven Hamiltonian systems [25-27], in which 
the time coordinate plays the role of an extra dimension, allowing for the realization 
of topological invariants with no equilibrium ground-state counterpart. In this work. 
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we focus on a different paradigm in which Hamiltonian unitary dynamics is replaced 
by specifically designed dissipative dynamics described by a quantum master equation. 
Such a scenario was originally proposed as a means for quantum state preparation 
and quantum computation [28, 29] and relies on the proper engineering of a coupling 
of the system to a suitable reservoir. In the context of cold atoms, such reservoir 
engineering may be seen as a natural extension of the more conventional Hamiltonian 
engineering, with similar advantages as compared to solid-state systems such as precise 
microscopic control and tunability. In previous works, we have shown how this concept 
can be utihzed to "cool" or drive ensembles of atomic fermions into topologically ordered 
states in one [30] and two [31] dimensions in a targeted way, starting from an arbitrary 
initial state described a density matrix. The analysis of the many-body properties of 
the phases and phase transitions arising in these examples has revealed similarities but 
also differences between the physics of topological ground states of Hamiltonians and 
topological steady states resulting from a purely dissipative evolution. 

In this work we put the results obtained in our two previous case studies into 
a broader theoretical perspective. We provide a framework for investigating non- 
equilibrium topological states that can be reached by means of engineered dissipation, 
developing a formalism and physical understanding that can also be used in situations 
where dissipation occurs as a perturbation. The natural object to study is the density 
matrix of the system, which docs not necessarily correspond to a pure state described 
by a wave function alone. In the present article we focus on quadratic master equations 
with the aim of classifying topological states described by density matrices in analogy 
to the Hamiltonian ground-state scenario. All information contained in the density 
matrix is then equivalently encoded in the covariance matrix gathering all static single- 
particle correlations. By identifying and exploiting the analogy between this object and 
a quadratic Hamiltonian in a "first-quantized" representation, we demonstrate how to 
classify topological phases in a non-equilibrium context where mixed states are allowed. 
Our analysis focuses on both bulk and edge properties. 

As compared to a Hamiltonian ground-state scenario, key differences arise from the 
fact that the dynamics — or the spectral properties of the system — and the properties 
of its "ground" (steady) state — or the static correlation properties — arc not as tightly 
related as in the Hamiltonian context. As far as the bulk is concerned, this crucial 
difference manifests itself in the fact that two independent spectral properties must be 
present to guarantee that the system is in a stable topological state: The first quantity 
that we identify is the dissipative gap, which corresponds to the slowest damping rate 
associated with modes belonging to the bulk of the system and is a direct counterpart 
of the excitation gap of a Hamiltonian spectrum. The second is the purity gap, which 
describes the purity of the mode belonging to the bulk which is most strongly mixed. 
Clearly, a purity gap is always present in the Hamiltonian context, since Hamiltonian 
ground states are by definition pure states. In our context, however, we argue that 
the system can undergo a topological phase transition if either (or both) of these two 
different gaps vanishes in a particular parameter regime. 
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The purity of the state plays a key role not only in the bulk, but also for the edge 
physics. In the Hamiltonian context, bulk-edge correspondence theorems describe a 
tight relation between the number of edge zero modes (i.e., modes that are decoupled 
from the Hamiltonian dynamics and thus do not evolve) found at the interface between 
two topologically distinct phases and the value of the topological invariant associated 
with each of the phases [13, 32-35]. We formulate a dissipative variant of such bulk-edge 
correspondence: Topological order ensures the existence, at the interface, of a fermionic 
subspace that is isolated from the bulk (with a dimension determined by the value of 
the topological invariant on both sides of the interface). However, in contrast to the 
Hamiltonian case, topology does not guarantee the decoupling of this subspace from 
the dynamics. As a result, the modes corresponding to this subspace can be either 
be zero-damping modes — i.e., modes that are decoupled from the dynamics similarly 
as in the Hamiltonian setting — or emerge as zero-purity modes — i.e., modes that are 
in a completely mixed state; in which no information can be stored. In the context 
of engineered dissipation, the simultaneous appearance of both zero-damping and zero- 
purity modes may give rise to intriguing physical effects, as we discuss in this work. 

The fact that actual physical implementations of model Hamiltonians need often 
be properly modelled as open systems due to particle losses or dephasing, e.g., has 
been recognized in a number of theoretical works focusing on the stability of the edge 
modes [36-39] or on the very definition of topological order in such circumstances [40]. 
We emphasize that our approach is fundamentally different here, since dissipation docs 
not occur as a perturbation but is rather harnessed as the main resource to generate 
the dynamics. 

Our paper is organized as follows. In section 2, we discuss the dissipative framework 
of interest: We introduce the concept of "dark states" in a many-body context, and 
explain the main ideas behind the physical implementation of a dissipative counterpart 
of Kitaev's quantum wire, thereby illustrating how to engineer more general dissipative 
evolutions giving rise to superfluid paired states. We also provide both a second- and 
a first-quantized formulation of quadratic dissipative dynamics, and discuss the key 
properties that are necessary to understand the bulk and edge physics of Gaussian 
states in terms of either the corresponding density matrix or the associated covariance 
matrix. In section 3, we construct a symmetry-based topological classification of 
driven-dissipative systems using the covariance matrix, and identify relevant topological 
invariants in one and two dimensions. In section 4, we then apply this framework to 
identify the classes D and BDI of Altland and Zirnbauer as the symmetry classes that can 
be dissipatively targeted under physical constraints related to "typical" implementation 
schemes. As is well known, the edge modes of systems belonging to these two classes 
are Majorana fermions, which explains the potential of dissipatively induced superfiuids 
to exhibit such modes. We also show that, in two dimensions, the quasi-locality of the 
dissipative operations acting on the system density matrix alone implies a vanishing 
Chern number. In section 5, we discuss the fate of the bulk-edge correspondence in the 
dissipative setting. We also show how to construct dissipative Majorana modes explicitly 
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in a translation-invariant setting, and examine the robustness of such modes in the 
presence of typical imperfections. Section 6 is devoted to the discussion of the physical 
role of the dissipative gap for adiabatic manipulations — in particular, for the braiding 
of dissipative Majorana modes — showing that dissipative Majorana modes exhibit non- 
Abclian exchange statistics just as their Hamiltonian counterparts. The remainder of the 
paper deals with illustrative examples of our general framework and results: In section 7, 
we analyze a "zigzag" dissipative quantum wire exhibiting topological phase transitions 
of the three possible types allowed by the closure of the dissipative and/or purity gaps. In 
section 8, we illustrate in a two-dimensional model a dissipative mechanism that makes 
it possible to obtain unpaired Majorana modes in a topological phase characterized by 
an even- valued integer topological invariant. 

2. Dissipative framework 

2.1. Quantum master equations for many-body systems 

The quantum master equation describing the time evolution of the reduced system 
density matrix p is given by 



The commutator term familiar from the Heisenberg equation describes the coherent 
dynamics generated by a system Hamiltonian H. The second part, often referred to 
as Liouville operator or Liouvillian, describes the dissipative dynamics resulting from 
the interaction of the system with an environment, or "bath". In particular, the set of 
Lindblad operators (or "quantum jump" operators) £i describe the coupling to that bath. 
The Liouville operator C has a characteristic Lindblad form: The anticommutator term 
describes dissipation and must be accompanied by fluctuations in order to conserve 
the "norm" Tr(p) of the system density matrix. The corresponding term, where the 
Lindblad operators act from both sides onto the density matrix, is called "recycling" or 
"quantum jump" term. Note that we have absorbed the rates associated with each 
dissipative process into the definition of the Lindblad operators, making them carry 
dimension of square root of energy. The rates are non-negative, so that the density 
matrix evolution is completely positive, i.e., the eigenvalues of p remain positive under 
the combined dynamics generated by H and C [41]. 

The quantum master equation provides an accurate description of a system-bath 
setting with a strong separation of scales. More precisely, there must be a fast energy 
scale in the bath (as compared to the system-bath coupling) that justifies to integrate out 
the bath in second-order time-dependent perturbation theory. If, in addition, the bath 
has a broad bandwidth, the combined Born-Markov and rotating-wave approximations 
arc appropriate, resulting in equation (1). Such a situation is generically realized in 
quantum optical few-level systems, e.g., for a laser-driven atom undergoing spontaneous 
emission. On the other hand, typical condensed matter systems do not display the 
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necessary scale separations to justify a microscopic description in terms of a master 
equation. In systems with engineered dissipation [42] as investigated in this paper, we 
are interested in scenarios that share features from both quantum optical systems — in 
that they are coupled to Markovian quantum baths — and condensed matter systems — in 
that they dispose of a continuum of spatial degrees of freedom on a lattice. Using the 
manipulation tools of quantum optics, the validity of a Markovian master equation can 
be ensured, giving rise to a well-defined microscopic dissipative many-body dynamics. 
A similar level of microscopic control as obtained in Hamiltonian engineering in a cold 
atom context can be expected for this "Liouvillian engineering", which therefore is a 
natural extension of the former to a more general non-equilibrium situation. In this 
context, dissipation does not occur as a perturbation, but rather as the dominant 
resource of the many-body dynamics. In particular, here we will consider the case 
where the Hamiltonian is absent H = 0. Such a scenario can be useful from a practical 
point of view — for cooling systems into desired states — but also gives rise to interesting 
new many-body physics. 

2.2. Dark states 

In the long-time limit, a quantum system governed by equation (1) approaches a 
stationary or steady state pf = p{t — >■ oo) which generically corresponds to a mixed 
state. An interesting situation appears if, instead, the many-body density matrix is 
driven towards a pure stationary state, pf = |'0)i3('0|r> [29,43]. In quantum optics, such 
pure states \iP)d that are obtained as a result of a dissipative evolution are called dark 
states. Mathematically, such dark states are zero modes of the Liouville operator. More 
precisely, a dark state is a zero mode shared by all Lindblad operators. 



The dark-state solution is the unique solution of the Liouville dynamics if (i) there exists 
exactly one normalized dark state \iP)d, and (ii) there are no stationary solutions other 
than this dark state [29, 44] . In the specific case of interacting Liouville operators 
(higher than quadratic in the creation and annihilation operators) [29, 44] or non- 
interacting Liouville operators (quadratic in the creation and annihilation operators), 
the fact that these conditions are satisfied can be proved rigorously. If present, the 
dynamics described by equation (1) for H = corresponds to a directed motion — in 
the Hilbert space of the system — into the "sink" provided by the dark state, which 
is reached independently of the initial density matrix. In recent years, a number of 
theoretical [28, 45] and experimental [46, 47] studies have focused on how to construct 
Liouville operators such that, in the long-time limit, a quantum system reaches a well- 
defined, pTirc many-body steady state or exhibits novel phase transitions resulting from 
the competition between coherent and dissipative dynamics [48-51]. In particular, in 
the context of atomic fermions, it has been shown how to engineer number- conserving 
dissipative dynamics that drives the system into a pure BCS-type paired state in the 
absence of conservative forces [52, 53]. The dissipative pairing mechanism forms a basis 
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Figure 1. (a) Unit cell for dissipation engineering. The lower potential wells 
correspond to the physical sites, whereas the upper site in between is an auxiliary site. 
Atoms in the lower sites are coherently coupled to the auxiliary site with opposite Rabi 
frequency ±fl. Decay back to the lower sites occurs via spontaneous emission, where 
energy is released into a surrounding reservoir (see text). If the coupling to the upper 
level is sufficiently far detuned (A ^ ^l), the latter can be integrated out, so that an 
effective dynamics for the lower sites is obtained, (b) In an optical lattice setup, this 
unit cell is repeated in a translation-invariant fashion multiple times in a natural way. 
The quantum wire corresponds to the lower sites — as anticipated above — which are 
populated by spinless (or spin-polarized) fermions in the cases discussed in this work. 
Rabi frequencies with alternating signs are realized by a commensurability condition 
on lattice and drive lasers (see text). Dissipation results from spontaneous phonon 
emission into a BEC reservoir (light grey). 



for the targeted cooling into states with non-trivial topological properties far from 
thermodynamic equilibrium. 



2.3. Physical realization 

Here we briefly sketch the implementation idea common to the dissipative models studied 
in this paper. The basic setting consists of an atomic ensemble confined in an optical 
lattice (the system), which is driven coherently and immersed into a larger reservoir 
consisting of a different atomic species and playing the role of the dissipative bath. In 
the case of interest here, the constituents of the system are fermions. In cold atomic 
gases, the associated spin is realized in terms of hyperfine states, and thus both the cases 
of spinless and spinful fermions can be meaningfully considered. The bath constituents 
are chosen as bosonic atoms, so that the conservation of fermion parity in the system is 
guaranteed. 

The working of the driven-dissipative mechanism is best illustrated by the unit cell 
A-configuration displayed in figure 1. The complete driven-dissipative process consists 
of two steps: The first step is a coherent excitation from the system (lower sites) to the 
auxiliary site in between. In the example of figure 1, we quasi- locally excite fermions 
on the system sites i and i + 1 (with annihilation and creation operators and aj 
corresponding to each site) into an antisymmetric superposition ~ — Oj+i, which can 



10 



be controlled by a commensurability condition of the driving laser to the standing wave 
laser generating the optical superlattice (see reference [42] for details): If the driving laser 
has twice the wavelength of the lattice laser, there is phase shift of vr in the effective Rabi 
frequency from one site to the next, leading to a relative minus sign. The auxiliary level 
is unstable if coupled to the reservoir. In this case, spontaneous phonon emission into 
the surrounding bath can occur, thereby giving rise to the second, dissipative step. The 
atoms are "brought back" to the lower sites in a quasi-local way ~ a]^ + al_^i] since this 
process is isotropic, there is now a relative plus sign. If this driven-dissipative process 
is generated using a drive laser that is far detuned (Q/A <^ 1) from the auxiliary site 
resonance frequency, the auxiliary site can be integrated out and we obtain a Lindblad 
operator of the form 

£j = CJA for alH, Cj = aJ + aJ+i, A = - a^+i. (3) 

A few remarks are in order: (i) For a driving laser superimposed over the extent of 
the whole optical superlattice, we obtain translation-invariant Lindblad operators as 
depicted in figure 1(b), up to system boundaries which are not shown, (ii) The quasi- 
locality of the operators is controlled by the Wannier function overlap between the 
onsite wave functions involved in the combined excitation and de-excitation processes, 
(iii) The Lindblad operators that can be realized in this setting have a generic form 
£. = cjAi, where Cj{Ai) is a linear translation-invariant superposition of creation 
(annihilation) operators with generic properties: The excitation part (A^) can be 
controlled to high accuracy — involving, in particular, the control of the relative phases 
in the superposition — since it is generated by a coherent laser beam. In two dimensions, 
for example, this allows to imprint angular momentum by choosing relative laser phases 
in different primitive directions of the lattice. On the other hand, the de-excitation or 
decay part (C|) results from spontaneous emission and is therefore unavoidably isotropic 
(or s-wave symmetric), (iv) The system particle number is conserved in our dissipative 
framework. This is reflected in the property [ii, iV] = for all Lindblad operators, 
where N = X]?'"'* ^^e total particle number operator, hi = a\ai. This exact microscopic 
property of the system, which implies an exact conservation of parity, is of importance 
to the discussion of the possible imperfections that may occur in the dissipative setting 
after performing approximations. Physically, this property originates from the fact 
that typical interactions in cold atomic systems are local density-density interactions. 
In particular, the system and bath constituents will interact via such coupling. On 
an even more elementary level, the fact that the bath is bosonic provides a further 
reason for fermion parity conservation. This aspect is crucially different from solid- 
state implementations which are not perfectly closed systems: There the environment is 
typically fermionic, which facilitates system-bath exchange processes affecting the parity 
of the system, (v) While no particle number exchange is possible between the system 
and the reservoir, energy can be exchanged. This enables the outflow of entropy from 
the system into the (inflnite) reservoir, and the targeted cooling into pure many-body 
states. A crucial prerequisite for the entropy removal from the system is the coherent 
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driving of the system, (vi) The fast energy scale ensuring the vahdity of the Born- 
Markov and rotating-wavc approximations underlying our construction is provided by 
the band separation between the system and the auxiliary sites, which is the largest 
energy scale in the problem, (vii) The creation and annihilation part of the Lindblad 
operators is respectively symmetric and antisymmetric under the exchange i <H- i + 1. 
Such property is important for reaching pure stationary states, as will become clear 
below. 

2.4- From interacting Liouvillians to quadratic master equations 

The physical implementation discussed above realizes a number-conserving microscopic 
dynamics, with the key advantage of conserving fermion parity as a consequence. The 
dynamics generated by the corresponding bilinear Lindblad operators is described by 
an interacting (quartic) Liouvillian. The Lindblad operators are constructed in such a 
way that the stationary state is a unique dark state for a given even particle number 
2iV, characterized by a BCS pair wave function with fixed particle number |BCS, N) 
that satisfies £j|BCS, N) = for all i. The Liouville operator ensuring this property 
thus represents a parent Liouville operator for a given fixed number BCS pair wave 
function (see Appendix A for more details). Starting from the exact knowledge of the 
fixed-number dark-state wave function, we can switch in the thermodynamic limit from 
a fixed-number to a fixed-phase ensemble. In particular, the long-time evolution of the 
interacting master equation can be linearized based on this knowledge. The calculation 
presented in Appendix A can be summarized as 



That is, the product of creation and annihilation parts in the quadratic Lindblad 
operators transforms into a sum, giving rise to linear Lindblad operators. This relies 
on the property that C|(Aj) is symmetric (antisymmetric). It provides a dynamical 
connection between the fixed- number and fixed-phase settings at the level of the equation 
of motion. The long-time dynamics is universal, in the sense that it is independent of 
the initial state. 

We note that a — re'^ is a complex number in the above equation. Its phase is 
not fixed by the dynamics, but rather reflects spontaneous symmetry breaking in the 
dissipative setting of interest. The modulus r, on the other hand, is determined by the 
average particle number in the system (see Appendix A). In particular, for half filling 
and in the example of equation (3), we find r — 1; such that for ^ = 0, without loss of 
generality. 



We recognize the quasi-local Bogoliubov quasiparticle operators of Kitaev's Hamiltonian 
quantum wire [15] (at half filling and with equal pairing and hopping amplitudes, up to 
normalization), emerging here naturally in the long-time limit of a dissipative dynamics. 
The ground-state condition of the Hamiltonian quantum wire, Lj|BCS, = for all 




(4) 




(5) 
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i, is now interpreted as the dark-state condition of the hnearized dissipative evolution. 
The corresponding wave function now has a fixed phase 9 instead of a fixed number. 
Since the Lindblad operators form a complete Dirac algebra, {Li^L}^] ~ 6ij, {Li, Lj} = 
{L], lJ} = for an infinite system with no boundaries, the uniqueness of the dark-state 
solution is obvious. 

The quadratic dynamics obtained in the long-time limit makes the systems 
considered in this work amenable to a treatment analogous to the discussion of 
quadratic Hamiltonians when examining their topological properties. This dynamics 
was obtained by giving up exact particle number conservation, which is justified in 
the thermodynamic limit. The absence of exact particle number conservation thus 
emerges similarly as in the Hamiltonian scenario. There is, on a formal level, however, 
a crucial difference between dissipative and Hamiltonian settings. While a quadratic 
number non-conserving BCS Ilamiltonian still conserves parity, formally such a property 
is not present in a quadratic Liouville evolution (for example, single fermions can be 
ejected into the environment, giving rise to quasiparticle poisoning [38,39]). However, 
remembering that the microscopic dissipative processes do conserve particle number 
and thus parity exactly, we can rule out parity non-conserving processes as potential 
imperfections arising in our scenario. The number non-conserving nature of the system 
is introduced "within the system only", and there is no exchange of particles with the 
reservoir. Physically, the number non-conserving processes describe pairwise creation 
and annihilation out of or into the mean field provided by the system itself. 

2.5. Gaussian master equations 

Having discussed how quadratic fermionic master equations naturally emerge in the long- 
time limit of interacting Liouville operators, we now summarize some general properties 
of such master equations. We do this in both a second- and a first-quantized formulation, 
working with operators or matrices, respectively, as familiar from the Hamiltonian 
setting [54]. For this discussion, it is useful to work in the real (or Majorana) basis 
of fermionic quadrature component operators. For a system with N sites, 2N real 
fermionic modes are introduced according to 

C2n-1 = i (On - 4) , C^n = 4, {c„. Cm} = 2Sn,m, 4 = C„. (6) 

The fact that the master equation is quadratic in the fermion operators implies solutions 
in terms of Gaussian density operators. In the second-quantized formulation, this can 
be written as 

Pit) ^ exp (^-'-c'G{t)cy (7) 

where = (ci, C2n) is a column vector defined from the 2N Majorana operators and 
G is a real antisymmetric matrix (so that iG is Hermitian). Formally, p thus has the 
form of a canonical density matrix pc ~ e~^^ for a quadratic Hamiltonian. 

Instead of working in second quantization, we can move to a first-quantized 
formulation. As in the Hamiltonian scenario, the latter allows us to discuss symmetry 
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and topological classifications in a more direct way. The key object here is the covariance 
or (equal-time) correlation matrix collecting the second moments, which is the only 
information contained in a Gaussian density operator. It is defined as 

We now look for an equation of motion for this object [49, 55]. A straightforward way 
to derive such an equation is via the adjoint equation to a given master equation for 
the density operator [56], describing the evolution of an operator in the Heisenberg 
picture. For the correlation operator r„^ = | Sj^=i ^^j^j'k^^k with real antisymmetric 
^jT = {hJm,k - Sj,mSn,k): this reads 

^ 2N 
dAm = Yl 4^nmLi - -{L\U, f„^} = -i C, ({M, 6?""^}),.fe Cfc. (9) 

i j,k=l 

Here we have written the linear quantum jump operators as Li = Ylk=i^k'^k, lI = 
X]fc=i ^fcCfe and introduced the matrix 

N 

Mjk^J2^ril = ^{X,k-'^Y,k), (10) 
1=1 

where X — and Y — —Y^ are real symmetric and antisymmetric matrices, 
respectively. Furthermore, X by construction is positive semidefinite. Taking the 
expectation value of equation (9), we readily find the fluctuation-dissipation equation 
describing the evolution of the real antisymmetric correlation matrix F, 

dtr^-{x,r} + Y, (11) 

where we have suppressed the matrix indices. Denoting the steady-state correlation 
matrix as F, which satisfies the equation {X, F} = Y, we can give a clear physical 
meaning to the matrices X and Y. Writing F = F -|- 5F, the approach to the steady 
state is governed by dtST — —{X,Sr}; i.e., X alone governs the damping dynamics 
towards that steady state. The matrix Y describes fluctuations, which come along with 
dissipation in a probability preserving (9tTr(p(i)) = 0) quantum mechanical evolution. 
The steady state F depends on both X and Y. 

Finally, we remark that the correlation matrix is related to the density operator 
equation (7) by 



^nm(t) = ^Tr (p(t)[cn, cj) = i tauh 



■ G{t) 



(12) 



2 

We may compare this to a Gaussian Hamiltonian equilibrium setting: Introducing the 
first-quantized, real and antisymmetric Hamiltonian matrix h in the Majorana basis via 
H — ^ "^ij hijCiCj, we have at an arbitrary temperature T — 1/P 

'./3h' 



Fi!S^ =itanh 



1 

2 



(13) 



which reduces to F*^®*^) = isgn(i/i) at T = 0. 
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2. 6. Purity and purity gap 

The purity of a Gaussian state defined by a particular correlation matrix F can be 
revealed by examining the spectrum of the Hermitian positive semidefinite matrix (iF)^, 
which we refer to as the purity spectrum. Pure Gaussian states are characterized by a 
"flat" purity spectrum with eigenvalues all equal to 1, whereas mixed Gaussian states 
exhibit eigenvalues smaller than 1, each zero eigenvalue indicating the existence of a 
completely mixed subspace. Intuition regarding the purity spectrum can be gained 
by constructing a fictitious quadratic Hamiltonian Hr from the Hermitian matrix iF, 
namely, 



where the q are the 2N Majorana basis operators introduced in the previous section. 
Since F is a real antisymmetric matrix, the spectrum of this Hamiltonian consists of real 
eigenvalues ±e„ {n — 1,2, . . . , N). Most importantly, the positive part of the spectrum 
of Hr is the purity spectrum of the Gaussian state represented by the correlation matrix 
F (up to a square root). Exploiting this analogy further, we will identify the eigenvectors 
of (iF)^ as "eigenmodes" or "modes" (of the fictitious Hamiltonian i^r)- In particular, 
we will refer to modes of Hr associated with zero eigenvalues as zero-purity modes and 
to the spectral gap of the latter as the purity gap. Such modes are defined in the 
mode space M. consisting of operators of the form 7 = v'^c with v G M?^ . Modes 
corresponding to a unit vector v will be referred to as "Majorana" modes since they 
satisfy 7^ = 7 and 7^ = 1. Moreover, we will distinguish two types of zero-purity modes: 
intrinsic ones, which are determined by the dissipative dynamics, and eoctrinsic ones, 
which result from mixed initial conditions (and thus disappear when starting from pure 
initial states). Prom a topological perspective, Majorana zero- purity modes that have a 
topological origin will be referred to as genuine, as opposed to spurious ones. 

The purity of the steady state is determined by the dissipative dynamics and, if the 
steady state is not unique, by the purity of the initial state (i.e., by the initial conditions). 
In the case of interest in this work where the dissipative dynamics is quadratic, one can 
show (we refer to our previous work [31] for an explicit proof) that there exist initial 
conditions leading to a pure steady state whenever the corresponding Lindblad operators 
Li form a set of anticommuting operators, i.e., whenever {Lj, Lj} = for all i, j *. In the 
matrix representation defined in the previous section, one can then establish a one-to-one 
correspondence between the matrices X and Y encoding the dynamics. Intuitively, this 
can be understood by examining the steady-state equation {X, T} — Y (where F now 
denotes the steady-state correlation matrix): if the steady state is pure, the spectrum 
of the associated correlation matrix F (i.e., the purity spectrum) essentially contains no 

*Note that Lindblad operators satisfying this condition generate the same (exterior) algebra as 
fermionic annihilation operators. They need not be fermionic annihilation operators, however. The 
anticommutation relation {Li, Lj} = (for all is a necessary and sufficient condition to ensure the 
existence of a pure state \tp) such that Li\ip) = for all i, which is all that we need. 




(14) 
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information, since all of its eigenvalues are equal to 1. The information contained in X 
must therefore be exactly the same as that contained in Y, otherwise the steady-state 
equation would not be satisfied. In other words, X and Y both contain full information 
about the dissipative dynamics when the steady state is pure. In that case, one can 
construct yet another useful object encoding all information about the dynamics: the 
so-called parent Hamiltonian i^parent naturally associated with the Hermitian matrix iy, 
defined as 



Clearly, the spectrum of -ffparent is directly related to that of Y and therefore to that 
of X as well for dissipative dynamics whose steady state is pure. Remembering the 
definition of the matrix Y in terms of the Lindblad operators, one can rewrite the 
parent Hamiltonian in the equivalent form 



This shows that pure steady states {ip), which are "dark states" satisfying Li\ip) = for 
all i (see equation (2)), can equivalently be seen as ground states of //parent- As opposed 
to the purely fictitious Hamiltonian ifp that we have constructed above to quantitatively 
assess the purity of an arbitrary Gaussian state, independently of any dynamics, the 
parent Hamiltonian therefore describes features associated with the actual (dissipative) 
dynamics of the system — as expected from its definition from the matrix Y. In fact, 
we will argue in the next section that the spectrum of i/parent encodes all information 
regarding pure steady states. We emphasize, however, that the parent Hamiltonian does 
not play such a prominent role in the more general case where the steady state of the 
dissipative dynamics is mixed (even when starting from pure initial states). 

As demonstrated in our previous work [31], the above discussion can be formalized 
and summarized as the following equivalent statements: 

(i) The steady state is pure (at least for pure initial states); 

(ii) {Li, Lj} = for all i, j (i.e., the Lindblad operators form a set of anticommuting 

operators); 

(iii) [X, Y] = and X"^ = — \Y'^ (in particular, the spectra of X and Y are directly 



(iv) The dissipative dynamics can be fully described using the parent Hamiltonian 

//parent ^i^i' 

This last point will be clarified in the next section. 

2. 7. Dissipative gap and Majorana zero-damping modes 

In the case where the dissipative dynamics is quadratic, the dynamical approach to 
the steady state is governed by the associated matrix X, as mentioned in the previous 
section. This matrix, which is by construction real, symmetric and positive semidefinite. 




(15) 




(16) 
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can be spectrally decomposed in the form X = Yl'j=i'^j{'^j ® ^J) with eigenvalues 
Kj > and associated eigenvectors Vj e M^^. The eigenvectors of X define "modes" 
in the mode space defined in the previous section. Assuming that they are normalized 
to unity, each eigenvector Vj can be identified with a corresponding Majorana mode 
jj = vjc. Physically, the eigenvalues of X then correspond to particular damping rates 
associated with particular Majorana modes. Accordingly, we will refer to the spectrum 
of X as the damping spectrum and to Majorana modes 7^ associated with a vanishing 
damping rate as Majorana zero-damping modes. One can show (for an explicit proof, 
see our previous work [31]) that a Majorana mode 7 = v-^c is a zero-damping mode 
whenever the following equivalent conditions are satisfied: 

(i) Xv = 0; 

(ii) {Li, 7} = for all i; 

(iii) if V ~ for all i (Ij being the vector corresponding to the Lindblad operator Lj 
in mode space, i.e., Lj = Ifc). 

We remark, however, that Majorana zero-damping modes satisfying the above conditions 
need not be spatially locahzed. The topological nature of the system will play a crucial 
role in ensuring such localization, as we will demonstrate in section 5.2 below. In general, 
we will distinguish genuine Majorana zero-damping modes that have a topological 
origin from spurious ones which a mere artefacts of the dissipative dynamics. A simple 
example of spurious modes is provided by a lattice site on which no dissipative dynamics 
takes place. This gives rise to two Majorana zero-damping modes decoupled from the 
dynamics, which obviously do not have a topological origin. 

The damping spectrum describes the dynamical separation (i.e., in time) between 
particular modes in a similar way as the spectrum of a Hamiltonian determines the 
energy separation between specific modes. Pushing the analogy further, one can see 
that the existence of a dissipative gap (or "damping gap") in the damping spectrum 
leads to the dynamical isolation of bulk and edge modes (through the quantum Zeno 
effect [57]), thereby providing a dissipative counterpart to the gap protection arising in 
the Hamiltonian context. Majorana zero-damping modes form a so-called decoherence- 
free subspace [58] unaffected by the dissipative dynamics and, in the presence of a 
finite dissipative gap, completely isolated from the rest of the system. The dissipative 
counterpart of a topological Hamiltonian ground-state degeneracy is then provided by 
the existence of a non-local decoherence-free subspace associated with zero-damping 
Majorana modes. 

We finally clarify the role of the parent Hamiltonian defined in the previous section 
in light of the considerations above. When the steady state of the dissipative dynamics is 
pure, the spectrum of X (i.e., the damping spectrum) directly maps to the spectrum of 
Y, which in turn trivially maps to that of ifparent- The parent Hamiltonian thus contains 
all information about the dissipative dynamics in that case. When the steady state 
is mixed (independently of the initial state), however, the tight relationship between 
X and Y (or i^parent) breaks down and the parent Hamiltonian becomes insufficient 
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to describe the dynamics. In this more general case, one can show (we refer to our 
previous work [31]) that a zero mode of i^parcnt (or, equivalently, of the matrix Y) does 
not necessarily correspond to a zero mode of X (i.e., to a zero-damping mode), although 
the converse is always true. In other words, zero modes of the parent Hamiltonian need 
not be Majorana zero-damping modes of the dissipative dynamics. In fact, any zero 
mode of i^parent which does not coincide with a zero mode of X gives rise, in steady 
state, to an intrinsic zero-purity mode. This crucial phenomenology will be key to 
understanding the non-equilibrium topological effects that will be exemplified below. 

To conclude this section, we remark that Majorana zero-damping modes do not 
benefit from the protection mechanism featured by i^parent due to the antisymmetry of 
the matrix Y. While the antisymmetry of Y forces i/parcnt to have an even number of 
Majorana zero modes, such that spatially isolated modes cannot be afi^cctcd by local 
perturbations, the fact that X is symmetric does not lead to such restriction. Although 
this is potentially harmful for Majorana zero-damping modes, we will demonstrate in 
the sections below that this can also lead to intriguing physics with no Hamiltonian 
ground state counterpart. 

3. Topological properties of the bulk 

In this section, wc focus on the topological properties of the bulk of drivcn-dissipative 
fermionic systems with Gaussian steady states. In particular, wc identify the correlation 
matrix, which fully describes such states, as a fictitious first-quantized Hamiltonian 
and use the latter to construct a topological classification in complete analogy to the 
conventional Hamiltonian scenario. 

The topological classification of gapped states of non-interacting fermions can 
be achieved on the basis of symmetry properties of the corresponding Hamiltonian 
under time-reversal, charge-conjugation (or particle-hole) and chiral (or sublattice) 
transformations, as was proposed by Schnydcr et al [59] following the classification 
of random matrices developed by Altland and Zirnbauer [54]. Ten symmetry classes 
were proved to be sufficient for an exhaustive classification of topological phases in 
any spatial dimension, and an alternative approach was later introduced by Kitaev in 
the powerful framework of topological K-theory [60-63] (see reference [64], e.g., for a 
self-contained review) . We argue below that all classification schemes developed in the 
Hamiltonian setting can be automatically applied to classify the Gaussian steady states 
of a dissipative dynamics. We do not provide an exhaustive classification, however, 
since this can be done straightforwardly based on the references cited above. Instead, 
we construct an explicit topological classification for two symmetry classes of particular 
interest for this work, namely, for dissipative systems belonging to the symmetry classes 
D and BDI. 
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3.1. Steady-state symmetries 

We first study how symmetries of the Lindblad operators translate into symmetries of 
the correlation matrix. To this end, we consider a Gaussian dissipative dynamics with 
unique steady state, i.e., the corresponding correlation matrix F is a unique solution of 

{x,r} = y, (17) 

with matrices X and Y defined as in section 2.5 (see equation (10), in particular). We 
assume that the Lindblad operators Lj are invariant, up to a phase factor, under some 
symmetry group G: 

g-^LiQ^e'^'Li, (18) 

where g E G. In order to preserve the linearity of the Lindblad operators in the 
fermionic operators, G must act linearly on the 2N Majorana operators Cj introduced 
in section 2.5 above. These operators form an orthonormal basis of the mode space 
J\A = of operators A = a^c with a e R^-^ with respect to the inner product 
{A, B) = [1/2){A, B} [65]. Clearly, any symmetry g E G (unitary or antiunitary) 
must act linearly on M and transform the Majorana basis defined by the operators cj 
into another Majorana basis. In mode space, a symmetry g & G must therefore be 
represented by a real orthogonal matrix Sg G 0{2N), 

g-'cg = SgC. (19) 

Note that this formula allows to analyze the symmetry properties of the state also in 
the more general case where the dynamics is governed by both a Liouvillian and a 
Hamiltonian, sec equation (24) below. 

The Lindblad operators are defined in the Nambu space J\f ^ C^^ of operators 
L — Fc with I e C^^, which can be viewed as a complexification ol M.. In Nambu 
space, the relevant representation of e G is given by Sg if g is unitary and SgK 
{— KSg) if g is antiunitary. Here K is the complex conjugation operator defined such 
that Ki = —iK. The Lindblad operators Lj = If c are therefore invariant (up to a phase 
factor) under the symmetry g & G if and only if 

li = e'^f^Sgli if g is unitary, (20) 

Ij = e"^'SgKli if g is antiunitary, (21) 

where 0i e [0,27r) (note that the phase factors e"^' do not affect the form of the 
LiouviUian). Using equation (10), we then find that the matrices X and Y encoding the 
dissipative dynamics have the properties 

X = SgXSj, Y = ±SgYSj, (22) 

where the positive and negative signs corresponds to the cases where g is unitary or 
antiunitary, respectively. The steady-state equation can then be written as 

{X, ±SjTSg} = Y (23) 

and, since we have assumed that the steady state is unique, we obtain 

r = ±sJrSg. (24) 
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This shows that the symmetries of the Lindblad operators are naturally reflected in 
symmetries of the steady-state correlation matrix T *. In turn, this implies that 
the matrix T can be used to construct a symmetry-based topological classification of 
Gaussian steady states. Of particular importance for this purpose are the two discrete 
symmetries corresponding to particle-hole (PHS) and time-reversal (TRS) symmetry, 
respectively. The former corresponds to the sign in equation (24) and is trivially 
satisfied as = 1 in our case, as we argue below, while the latter corresponds to the 
"-" sign and depends on the specific form of the Lindblad operators. 

We remark that chiral symmetry (defined as the combination of PHS and TRS [33]) 
is automatically satisfied whenever TRS is present, since the system always has PHS, 
by construction. In that case, there exists matrices Sc and St corresponding to PHS 
and TRS, respectively, such that equation (24) is satisfied (with a "+" sign for Sc and 
a "-" sign for St), and one can easily verify that the combination of PHS and TRS (i.e., 
chiral symmetry) leads to {ScS^, T} = {StS^, P} = 0. 

3.2. Topological classification and topological invariants 

Let us consider an arbitrary Gaussian steady state p, fully characterized by its 
correlation matrix Pjj = (i/2) Tr(p[Q, Cj]). Using that the matrix P is real and 
antisymmetric, we construct a corresponding fictitious free-fermion Hamiltonian as in 
section 2.6 where the purity spectrum was defined (not to be confused with the parent 
Hamiltonian introduced in the same section). 



thereby establishing a one-to-one correspondence between the set of Gaussian steady 
states and the set of free-fermion Hamiltonians. It is clear that the symmetries of 
the dissipative system — embedded in P — are the same as that of the Hamiltonian 
system defined by Hr, since P can be viewed as the "first-quantized" Hamiltonian 
corresponding to the "second-quantized" Hamiltonian Hi. The problem of classifying 
steady states according to topological properties is therefore equivalent to that of 
classifying Hamiltonian systems of non-interacting fermions, which is the main message 
of this section. Consequently, both the symmetry-based classification of references [54, 
59] and the K-theory approach of references [60-62] can be directly apphed in the 
dissipative framework. Note that the Hamiltonian Hr always takes a Bogoliubov-de 
Gennes form when expressed in terms of the original fermionic operators a^, a|, and is 
therefore automatically particle-hole symmetric. 

The topological classification crucially relies on the existence of a bulk spectral 
gap and is essentially based on the mathematical concept of homotopy equivalence, 

*Clcarly, using equations (8) and (19), the transformation of the correlation matrix is F — > zbS'jFS'g 
for unitary or antiunitary symmetries, respectively, irrespective of the dynamics — which in particular 
may involve both Hamiltonian and Liouvillian parts. In contrast, here we focus on how the invariance 
of the Lindblad operators (under some symmetry) translates as that of the correlation matrix. 




(25) 
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or equivalence under continuous deformations *. More specifically, two gapped 
Hamiltonians and H^' are considered as "topologically equivalent" if they can 
be continuously deformed into each other without closing the gap. The dissipative 
counterpart of this equivalence is provided by the mapping of equation (25): two 
Gaussian steady states corresponding to correlation matrices T and T' will be considered 
as topologically equivalent and referred to as belonging to the same topological phase if 
and only if they can be continuously deformed into each other without closing the bulk 
purity gap f- The existence of a bulk purity gap is therefore the key ingredient required 
to define topological order in the dissipative setting. 

The mapping defined by equation (25) and the results of references [59-62] provide 
us, in principle, with a general topological classification of all possible (purity) gapped 
Gaussian steady states according to symmetries and to the spatial dimension of the 
system. We now sketch this construction focusing on the symmetry classes that are 
most relevant for the dissipative systems considered in this work. 

An arbitrary 2N x 2N steady-state correlation matrix F {N being the number of 
fermionic modes, or the number of sites for systems of spinless fermions defined on a 
lattice) can be brought to a block diagonal form 



where Q is an orthogonal matrix and ±e„ are the real eigenvalues forming the spectrum 
of the Hermitian matrix iF. The purity spectrum is defined by the spectrum of the 
real symmetric matrix (iF)^, which is doubly degenerate with positive eigenvalues e^. 
Assuming that it is gapped, such that |e„| > for all n, the matrix F can be continuously 
deformed into a topologically equivalent matrix f with a "flat" purity spectrum 



Since (iF)^ = 1, this "spectrally flattened" correlation matrix defines a pure Gaussian 
state which is topologically equivalent to the not necessarily pure original steady state of 
interest. The matrix F allows us to construct an orthogonal spectral projection operator 
P (see, e.g., reference [13]) defined as 



which projects onto the A'"-dimensional subspace associated with eigenvectors of iF 
with negative eigenvalues %. This operator plays a crucial role in the topological 

*In the framework of K-theory, stable equivalence also plays a crucial role in the definition of 
"topological equivalence" (see reference [60]). 

fNote that the spectrum of -ffr is defined by that of the Hermitian matrix iF, and is therefore in 
one-to-one correspondence with the spectrum of (iF)^, i.e., with the purity spectrum. 

|From the point of view of the Flamiltonian H-p (see equation (25)), P projects onto the subspace 
of eigenstates of with negative energy. 




(26) 




(27) 



^(i-if). 
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classification of Gaussian steady states as well as in the construction of associated 
topological invariants, as we will demonstrate below. 

We first consider the case of spinless fermions on a d-dimensional lattice 
with periodic boundary conditions evolving under a translation-invariant dissipative 
dynamics. It will be convenient to label the local Majorana operators q as Ci^x, where 
i refers to a particular lattice site at position rj and A = 1,2 distinguishes the two 
local Majorana operators associated with the corresponding fermionic creation and 
annihilation operators Oj and a|, i.e., q^i — a] + Ui and Ci^2 — — Oi)- In momentum 
space, the steady-state correlation matrix F then takes the form of a 2 x 2 complex 
antihermitian matrix r(k) with components {X, /i — 1,2) 

(r(k))A, = ^ e''^-('-^-^')ra,, = Tr(p[ck,A, c_k,,]), (29) 

where N is the total number of lattice sites and Ck,A = G~^^'^'Ci^x- Since the matrix 

r(k) satisfies the condition r(k)* = r(— k), the spectrum of the Hermitian matrix ir(k) 
is given by eigenvalues ±e(k) with e(k) = e(— k) > 0. This allows us to introduce 
the "spectrally fiattened" momentum-space correlation matrix r(k) (the eigenvalues of 
ir(k) being ±1) and the associated spectral projection operator 

P(k) = i(l-if(k)) (30) 

projecting onto the subspacc C^(k) of the 2-dimcnsional complex vector space 
spanned by the (complex) eigenvectors of ir(k) associated with the negative eigenvalue 
— 1. The operators P(k) form a family of operators labeled by the wavevectors k 
belonging to the Brillouin zone, i.e., to the d-dimensional torus T'' = [— tt, tt]''. They 
define a fiber bundle ®ker<^C-(k) on the Brillouin zone manifold, with fibers C^(k) 
assigned to each point k 6 T^. The problem of classifying Gaussian steady states 
according to topology therefore reduces to that of classifying fiber bundles over a torus 
T*^ (the Brillouin zone). In the present case, the fibers are vector spaces and the fiber 
bundle is thus a vector bundle, for which a complete topological classification can be 
constructed using the approach of reference [59] or K-thcory [60-62] . Wc present below 
the results pertaining to the two types of systems that are most relevant to this work 
(see section 4), namely: 2D dissipative systems without TRS (symmetry class D of 
Altland and Zirnbauer) and ID dissipative systems with TRS (symmetry class BDI). 

The topological classification generally reduces to the identification of all possible 
homotopy classes of continuous maps k i— ?> P(k) from the Brillouin zone manifold to the 
space of spectral projection operators. Since the 2x2 complex matrix ir(k) is Hermitian 
and unitary, with (ir(k))^ — 1, the spectral projection operator can be expressed as 

P(k) = i(l + n(k)-cT), (31) 

where n(k) is a unit vector on the 2-sphcrc (S^) and cr is a vector of Pauli matrices. In 
the absence of additional symmetries, the spectral projection operator thus describes a 
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mapping from the Brillouin zone torus to the sphere *. 

In 2D dissipativc systems where the steady state belongs to the symmetry class 
D, the only symmetry that r(k) (or -P(k)) has is the PHS automatically embedded in 
Gaussian steady states. The homotopy classes of maps k — > n(k) then form a group 
which can be identified with the homotopy group 7T2{S'^) = Z f (which is non-trivial) and 
one can have non-trivial topological steady states and distinguish them via an integer 
topological invariant known as the Chern number, defined as 

dP dP dP dP 



1^20 = ^ y J dk^dky Tr 



dkrdL, Tr 



dkj^ dky 0ky dkr/; 

dr or or or 



167r7_^7_^ ^ y y \dkxdky dkydk^ 

=i/:/_>-»(-(t;S)' 

where in the last hue we have used the fact that r(k) = i(n(k) • cr) (see equation (31)). 

In ID dissipative systems where PHS is the only symmetry, the Brillouin zone 
corresponds to a circle and the homotopy classes of maps k i— )■ n(A:) form a homotopy 
group 7ri(S'^) which is trivial. This means that the vector n(A;) cannot "twist" in 
a way that cannot be continuously untwisted as we move along the Brillouin zone 
circle, such that all steady states are necessarily topologically trivial. The situation can 
change in the presence of additional symmetries, however, since additional constraints 
on P(k) (or, equivalently, on f (k) or n(k)) can potentially restrict the set of allowed 
continuous deformations. If the ID dissipative system has TRS, then it also has chiral 
symmetry and, hence, belongs to the symmetry class BDI. In that case, as argued below 
equation (24), one can find a unitary matrix E such that = 1 and 

{E,f(k)} = (33) 

for all k. (In our case, TRS takes the form Cj,i — > Cj,i, Ci,2 — >■ — Ci,2, such that E = a^.) 
The matrix E has eigenvalues ±1 and can be expressed as E = a-<T, where a is a real unit 
vector that does not depend on k. The above condition then takes the form a • n(k) = 
for all k. In other words, chiral symmetry restricts the vector n(k) from the sphere 
to a circle in the plane perpendicular to a. As a result, the relevant homotopy group 
becomes non-trivial, given by tti^S^) = Z. The corresponding topological invariant can 
be constructed by choosing a basis where 

with e(k) e [0,27r) and e'^C*) = %(k) + in^(k). The U(l) phase ^(k) therefore fully 
characterizes r(k) and contains all required information to construct a topological 

*P(k) can alternatively be seen as an element of U{2)/U{1) x f7(l) (7(1,2), G{m,n) being the 
set of all n-dimcnsional subspaccs in C™ (also known as a complex Grassmann manifold) [59]. 

t7r2(S'^) formally classifies maps from S"^ to S^, but since 7ri(S'^) is trivial, it equivalently classifies 
maps from to S"^. 
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invariant distinguishing between different topological classes of steady states. The 
relevant topological invariant, the winding number, takes the form 

2 r..^^.o 



— — I dfc a- nx — 

TT ./_^ V V dfc 



dfc 



where we have introduced projection operators P± = (1 ± S)/2 defined so that 

p,r(k)P_ = ( » ) , P_r(k)p, = ( _^ ° ) . (36) 

We remark that the above construction of topological classes and topological invariants 
crucially relies on translational symmetry. However, one expects bulk topological 
properties to be robust against weak dissipativc perturbations that preserve the purity 
gap as well as the symmetries of the system. This can be rigorously demonstrated 
in the K-theoretic framework of reference [60] and in the more general framework 
of noncommutative geometry developed in references [13,66], allowing to extend the 
above construction to disordered (and possibly finite) systems. Although the details are 
beyond the scope of this work, we emphasize that the existence of quantized (integer) 
topological invariants in such systems crucially relies on the quasi-local nature of the 
spectral projection operator. Here the quasi-local nature of the corresponding steady- 
state correlation matrix is ensured in the presence of a dissipative gap, which is a spectral 
property of the model in contrast to the purity gap, which relates to mode occupations. 
More specifically, the existence of a finite dissipative gap A implies an exponentially fast 
relaxation to the steady state (on a characteristic time scale r ~ 1/ A) and leads to the 
exponential decay of all spatial correlations, such that the spectral projection operator 
(see equation (28)) satisfies 

\PiX,j^ ^ cexp(-Q;|rj - r^l) (37) 
for some constants c, a > *. Note that the converse is also true [50]. The operator 
defined by {Pij)\^ = Pixjfj, (see reference [13]) is thus quasi-diagonal, i.e., there exists 
some constants c', c" > and a' > d {d being the spatial dimension of the system) such 
that 

\\Pij\\m < c'|r,-r,-|-"', ||P,,|| < c", (38) 

where ||.||hs is the Hilbert-Schmidt norm and ||.|| the usual operator norm. This 
"localization" condition satisfied by the spectral projection operator in the presence 
of a finite dissipative gap crucially allows for the definition of the Chern number and 
winding number topological invariants in disordered finite systems. Such a construction 
can be found in the appendix C of reference [13], for example. 

*Notc that the purity gap of F ensures that the exponential decay of all spatial correlations remains 
under the continuous "spectral flattening" transformation F — >■ F. 
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3.3. Phenomenology of non- equilibrium topological phase transitions 

The above discussion shows that the existence of two gaps is necessary to identify the 
topological properties of the bulk (Gaussian) steady state: (i) the purity gap, which 
allows to identify a particular topological class (e.g., Z) to which the steady state 
belongs, and (ii) the dissipative gap, which ensures that the steady-state correlations 
are quasi- local (as defined above), so that in particular the Chern and winding number 
topological invariants defined above (sec equations (32) and (35)) remain quantized 
in the presence of disorder. Since the dissipative and purity spectra are in general 
independent quantities — which can be seen from the fact that two independent matrices 
X and Y are required to describe the dissipative dynamics — ^these two gaps can close 
independently of each other. However, the closure of the dissipative gap alone gives 
rise to a critical behavior in steady state. As a consequence, non-equilibrium topological 
phase transitions can occur in three distinct ways: 

(i) Via the closure of the bulk dissipative gap only (criticahty) ; 

(ii) Via the closure of the bulk purity gap only (no criticality ) ; 

(iii) Via the closure of the both the bulk dissipative and purity gaps (criticality) . 

These three possibihties will be exemplified in exphcit models in section 7 below. Clearly, 
the same phenomenology must appear at interfaces between distinct non-equilibrium 
topological phases. We will investigate such physics in detail in section 5 below. 

4. Physical constraints 

Our discussion so far was very general, with the only assumption that the steady state 
of the dissipative dynamics is a Gaussian fermionic state. We now discuss restrictions 
that arise in "typical" experimental reahzations, as sketched in section 2.3 above. More 
specifically, we consider quasi-local dissipative processes in one- and two-dimensional 
lattice systems; assuming that the corresponding Lindblad operators act on a quasi- 
local subset of sites and have a spatially isotropic (or s-wave symmetric) creation part 
resulting from spontaneous emission processes. We additionally assume that the system 
is translation-invariant and that the steady state of the dissipative dynamics is unique 
and pure. Under this purity assumption, the steady state of the system can be viewed as 
the ground state of the parent Hamiltonian associated with the dissipative dynamics (see 
equation (15) and the discussion of section 2.6), and the existence of topological order 
can be assessed in the exact same way as in the Hamiltonian setting. In particular, 
as for the topological classification of section 3.2 above, one can rely on the general 
Hamiltonian classification of gapped topological phases of non-interacting fermions 
(according to symmetry and spatial dimension) developed in references [33,54,59,67]. 

Although we assume a pure steady state, we remark that the conclusions drawn 
below will also hold in the more general case where the steady state is not necessarily 
pure but has a gapped purity spectrum. We have argued in section 3.2 that an arbitrary 
state exhibiting a finite purity gap is topologically equivalent to a pure state exhibiting 
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a completely "flat" purity spectrum. To identify the possible classes of topological 
steady states that can be reached when taking account the "typical" physical constraints 
described above, one can therefore focus exclusively on pure states, without loss of 
generality. We emphasize that in that case the (automatically "flat" ) correlation matrix 
ir describing the pure steady state precisely corresponds to the parent Hamiltonian 
in its "first-quantized" form, namely, i^parent = ^^ij^ijCiCj = Hf, where Hf is the 
fictitious Hamiltonian that can always be constructed from the (Gaussian) steady state 
(see equations (25) and (28) and discussion thereof). 

Translational symmetry makes it most convenient to express the Lindblad operators 
in the momentum-space form = ttkOk + ^kfllk- order for our assumption of a pure 
and unique steady state to be satisfied, the Lindblad operators must form a complete 
set of anticommuting operators (i.e., {Lk,Lk/} = for all k, k'), which implies that 

UkVk = -li-k^^-k- (39) 

Translational symmetry additionally restricts the possible lattices for the system to 
Bravais lattices. In that case, the set of Lindblad operators is complete whenever there 
are as many Lindblad operators as lattice sites. 

As argued above, the pure steady state can equivalently be described as the ground 
state of the parent Hamiltonian i^parent = Sk-^k-^k (see section 2.6), which takes here 
the explicit form 

-f^parent = X](«k«k + ^^k«-k)(«kak + VkO^-k) 



= 2Z(|Mk|^ - \v-kf)alak + ^ |wk|^ 

k k 

+ IY. [("k^'k - «*_k^-k)aLaLk + h.c^ ■ (40) 

k 

Dropping the constant term \ it becomes 



-f^parent = ^ ^kCtlcik + ^ (^k^L^-k + ^k^-kOk) 



(41) 



k 

where we have defined 

Ck = kkl^ - k-k|^, Ak = u^Vk - lilk^-k = A_k. (42) 

Since the Lindblad operators are generally defined up to a phase, we can assume to 
be real, without loss of generality. Using equation (39), we thus find 

Ak = 2uivi,. (43) 

Our assumption that the creation part of the Lindblad operators is isotropic implies 
that Vk = V-k *• Since equation (39) must be satisfied, we conclude that the condition 

*In general, it can occur that = V-^ is only satisfied up to a phase e''^^''^ where 0k is a linear 
function of k, since the creation part of the Lindblad operators need not be isotropic with respect to 
a center of symmetry that corresponds to a lattice site. In that case, however, uu carries the opposite 
phase — since the annihilation part of the Lindblad operators is engineered with respect to the same 
origin (see section 5.2) — and our analysis remains the same, without loss of generality. 



26 



Uk — —U-k must be imposed. In practice, this can easily be achieved by tuning the 
phases of the driving lasers, as argued in our previous work [31], and will therefore be 
assumed in the following. We then have 

= C-k; Ck = M - (44) 



4.1. Symmetry classes 

Equations (41), (43), and (44) show that the parent Hamiltonian takes the generic form 
of a Bogoliubov-de Gennes (BdG) Hamiltonian and therefore exhibits, by construction, 
particle- hole (or charge-conjugation) symmetry (PHS). In the scheme of reference [33], 
further classification can be achieved by considering time-reversal symmetry (TRS). In 
the Nambu representation, the parent Hamiltonian can be written as 



^^parent ^ ^l^k*k with T^k 



2 
k 




(45) 



where is the "first-quantized" parent Hamiltonian and ^^j^, = (a|^, a_k) the Nambu 
field operator. Particle-hole and time-reversal symmetries are then present whenever 
there exists 2x2 unitary matrices Uc and Ut, respectively, such that the following 
conditions on the Hamiltonian are satisfied: 

Uhn*_kUc = -T^k, (PHS) (46) 
U\H*_kUT = +?^k- (TRS) (47) 

The symmetry class of the system is determined by whether or not such matrices can 
be identified and, if so, on the sign of UqUc and U^Ut (which can only take values 
±1, as argued in reference [33]). In the present case, the parent Hamiltonian describes 
a spin-polarized superfluid and exhibits — again, by construction — PHS with Uc = (Jx, 
such that UqUc = +1- It is clear from the form of //parent (see equations (41), (43), 
and (44)) that the system additionally exhibits TRS if and only if Al^ = ^k which, 
remembering equation (42), is satisfied provided that Ak is real up to a global phase *. 
In that case, equation (47) holds for Ut = I (such that U^Ut = +1) and the system 
exhibits chiral symmetry. 

In summary, two symmetry classes of steady states can be realized through "typical" 
quantum-reservoir engineering (see section 2.3) in driven-dissipative systems of spin- 
polarized fermions: if the Lindblad operators break TRS, the steady state belongs to 
the symmetry class D of Altland and Zirnbauer; otherwise it belongs to the symmetry 
class BDI. According to the topological classification of reference [33] , one must therefore 
consider either (i) 2D systems in the symmetry class D (with broken TRS) or (ii) ID 
systems in the symmetry class BDI (with TRS) in order to have the (a priori) possibility 
of reaching topologically non-trivial steady states. In these systems, distinct topological 
states are characterized by distinct values of an integer-valued (Z) topological invariant: 

*One must keep in mind that the gap function Ak is defined up to a global phase, since -ffparent is 
invariant under the global U{1) gauge transformation aj^ e'^^iL 
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the Chern and winding number, respectively. Most importantly, the zero-energy edge 
modes appearing in such systems in the Hamiltonian setting are Majorana fermions. 
This opens up the possibility to generate spatially isolated Majorana zero- damping 
modes in the dissipative setting of interest in this work. 



4-2. Chern number 

Let us now examine more closely the case of a 2D translation-invariant dissipative 
system with broken TRS, with a pure and unique steady state that accordingly belongs 
to the symmetry class D of Altland and Zirnbauer and is characterized by the Chern 
number topological invariant of equation (32). We show below that quasi-local Lindblad 
operators generating such dissipative dynamics do not allow to obtain phases with a non- 
zero Chern number, in contrast to the Hamiltonian setting. Although this result seems 
to be a no-go statement for Majorana zero modes in such systems, we will demonstrate in 
section 8 that unpaired Majorana zero-damping modes can exist in dissipative systems 
with vanishing Chern number provided that the steady state is not pure but mixed. 

As we have discussed above, a pure steady state can equivalently be viewed as the 
ground state of the parent Hamiltonian -ffparcnt associated with the dissipative dynamics, 
given by equation (45) (see also equations (41), (43) and (44)), and its bulk topological 
properties can be characterized by means of the spectral projector Pk = |(1 + nk ■ cr), 
where • cr — "Hk/A^k (Ak being a normalization factor defined such that ||nk|| = 1). 
The corresponding Chern number is determined by the vector nk (see equation (32)), 
which takes the explicit form 

/ 2^;kRe(^tk) \ 

2vklm{uk) , (48) 
1^ |tikP - |vkp J 

where A^k = v^Ck + l^kp = \/|Mkp + I't'kp *■ Note that A^k corresponds to the spectrum 
of i^parent and thus defines the dissipative gap (see section 2.7). In order for the vector 
life to be well defined, the functions i^k and must not vanish simultaneously anywhere 
in the Brillouin zone. Introducing (p^ — |</7k|e'^'' = 'fk/i^k and a "Fermi surface" 
— {k |(/9k| > (defined for some arbitrary e > 0) which generally f takes the 
form of a finite union of piecewise smooth, simple closed curves J-'^^x (i.e., J-^ — (Ja •^e,x), 
the Chern number can be expressed (see reference [31]) as a sum 

A 

*Wc assume that the system is infinite such that the components of nk are continuous functions 
of k in the Brillouin zone, in which case the Chern number given by equation (32) is a well-defined 
quantity. 

f More pathological Fermi surfaces can in principle be encountered, but the Chern number is a mere 
definition in that case. 



1 / Re(Ak) 
-lm(Ak) 



v 



1 

Ak 
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of winding numbers defined by 

W,,x = Vk^k -dk^^ (f {dkji^dh + dk^e^^dky). (50) 

Choosing e <^ 1, the above expressions tell us that the value of the Chern number can 
be inferred from the behavior of </?k = ^k/^k around the zeroes of — which do not 
coincide with zeroes of since A^k must be non-zero. Since the norm of (/^k diverges 
upon approaching such points, it is the phase winding of (/9k around these points that 
encodes all information about the Chern number. More specifically, each of the zeroes 
of Uk contributes to the Chern number by an integer value corresponding to the winding 
number of the phase ^k around the latter (in a conventional, e.g., counter-clockwise 
direction). Since the function can always be chosen as real (see discussion below 
equation (43)), ^k corresponds to the phase of (up to a sign that can be absorbed 
into the definition of the Chern number by a change of convention z/2d — ^ —i^2d) and 
the Chern number simply reduces to the sum of the phase windings of the function ui^ 
around its zeroes. 

In the Hamiltonian setting, the phase of -Uk is defined by the phase of the gap 
function Ak, whereas the norm of is defined by the norm of Ak and the value of 
the independent quantity ^k (the single-particle dispersion) [12], such that the phase 
and the modulus of are essentially independent. In our dissipative setting, however, 
Uk is a single complex function defined on the Brillouin zone and, therefore, the phase 
and the norm of are closely related. This crucially restricts the possible values of the 
Chern number: the function U]^ defines a smooth vector field (Re(Mk), Im(itk)) on a torus 
and, according to the Poincarc-Hopf theorem, the sum of its phase windings around its 
zeroes must vanish. We thus conclude that the only value of the Chern number that 
can be achieved in the above dissipative setting is zero. 

Note that the vanishing of the Chern number is a direct consequence of the 
smoothness of the vector field corresponding to u^, which in turn is due to the quasi-local 
nature of the Lindblad operators. The function u^ defines a smooth (C°°) vector field if 
and only if the norm of the coefficients uj^i — u{rj — Rj) encoding the annihilation part 
of the Lindblad operators (see equation (52)) decays faster than any power of l/|rj— Rj|. 
A power-law decay of the Lindblad operators would therefore be necessary to be able 
to engineer phases characterized by a non-zero Chern number in our dissipative setting. 

5. Dissipative edge physics 

Thus far we have concentrated on bulk properties of purely dissipative systems described 
by a quadratic Lindblad master equation. It is clear, however, that the existence of non- 
trivial topological bulk properties must somehow be reflected in the physics occurring 
at physical edges or in topological defects — which we both refer to as "edges" in the 
following. In the Hamiltonian setting, this intimate connection has been rigorously 
demonstrated and formalized — at least, for non-interacting systems — as a bulk- edge 
correspondence (or bulk-boundary correspondence) relating the topological nature of 
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the bulk to the number of gapless modes occurring at an edge [13, 32-35]. Of particular 
interest here is the bulk-edge correspondence pertaining to (i) ID systems characterized 
by a winding number topological invariant z/id (see equation (35)) and (ii) 2D systems 
characterized by a Chern number topological invariant z/2d (see equation (32)), which 
tells us that a number m = jz/*^^^ — of gapless edge modes must be present at the 
interface between two topological phases characterized by integer topological invariants 
(winding or Chern numbers, as appropriate) u^^^ and respectively. In general, such 
modes are robust (e.g., against disorder) as long as the Hamiltonian system remains 
in the same topological class. Since the purity gap can close in the dissipative setting 
of interest in this work, it is highly unclear whether similar universal signatures of 
bulk topological properties can be observed. We thus proceed, in the next sections, 
to investigate the edge physics that may arise in systems described by a quadratic 
dissipative dynamics. 

5.1. Dissipative bulk-edge correspondence 

Let us consider a generic dissipative dynamics characterized by a finite dissipative gap 
and described by Lindblad operators Li or, cquivalently, by matrices X and Y as defined 
in section 2.5. It is clear that the bulk-edge correspondence pertaining to Hamiltonian 
systems must also be satisfied when the steady state of the dissipative dynamics is 
pure, since such a state can equivalently be regarded as a ground state of the parent 
Hamiltonian ifparent = X^j-^l-^i — iX^ij^jCjCj (see section 2.6). In that case, the 
spectrum of i^parent defines the damping spectrum, and at least m — Iv^^^ — v^'^^] 
Majorana zero-damping modes must be present at the interface between two non- 
equilibrium topological phases characterized by Chern or winding number topological 
invariants u'^^^ and z/*^^\ In the more general situation where the steady state exhibits a 
finite purity gap but is not pure, the Hamiltonian bulk-edge correspondence obviously 
still applies to the parent Hamiltonian, such that i^parent still exhibits a minimum of 
m — — Majorana zero modes at the interface, but none of these modes is 
guaranteed to be Majorana zero-damping modes of the dissipative dynamics. As argued 
in section 2.7, Majorana zero mode of iJparent must either correspond to Majorana zero- 
damping modes or give rise to intrinsic Majorana zero-purity modes in steady state. The 
dissipative counterpart of the Hamiltonian bulk-edge correspondence can therefore be 
stated as follows: As long as the topological nature of the steady state remains the same 
in the bulk — which is true if the purity gap remains finite and if the symmetries of the 
Lindblad operators are preserved — a minimum total number of m = lu^^^ — u^'^^ = Az/ 
Majorana zero-damping and intrinsic zero-purity modes must be present at the interface 
between two non- equilibrium topological phases characterized by topological invariants 
i/^^^ and respectively. With obvious notations, this result can be summarized in the 
form 

(^damping) edge + ("^purity) edge ^ {^^)bulk- (51) 
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We remark that such a dissipativc bulk-edge correspondence is only guaranteed to hold 
in the presence of a finite dissipative gap, since the existence of such a gap is necessary 
to ensure the exact quantization of the topological invariants, as discussed in section 3.2. 
Most importantly, the inequality appearing in equation (51) originates from the fact that 
"spurious" Majorana zero modes may arise "accidentally" , as opposed to Majorana zero 
modes which we refer to as "genuine" that have topological origins. A strict equality is 
therefore only obtained when restricting to genuine Majorana zero-damping and zero- 
purity modes. 

Despite its important predictive power, the above result does not provide specific 
information regarding the exact number of genuine Majorana zero-damping modes that 
appears at an edge — as we will argue below, such knowledge must come from additional 
considerations unrelated to the topological nature of the system. One may therefore 
wonder whether this number can be a robust quantity and, if so, under what conditions. 
This will be the focus of the next sections: we will first investigate systems whose edge 
dissipative dynamics emerges in a very natural way by extending the bulk dissipative 
dynamics as close as possible to a physical edge. Since the topological nature of generic 
Gaussian states is only well-defined in the presence of a finite purity gap and since all 
Gaussian states featuring such a gap can be continuously deformed to pure states (see 
section 3.2), we will focus on dissipative dynamics driving the system into a pure steady 
state. In that case, the purity of the system obviously forces the number {m^^T-ity) edge 
of Majorana zero-purity modes to vanish in equation (51), and one can try to construct 
Majorana zero-damping modes explicitly. Having examined this simple situation where 
the number (mdampmg)edge of genuine Majorana zero-damping modes is solely determined 
by the topological nature of the bulk, we will then extend our discussion to the more 
general case where the dissipative dynamics can affect the purity of the system while 
maintaining a finite purity gap in the bulk, and conclude by clarifying the role of 
topology in the dissipative setting central to this work. 

5.2. General form of Majorana zero-damping modes for translation-invariant 
dissipative processes 

We have introduced above a dissipative bulk-edge correspondence which crucially rehes 
on the existence of two gaps: the dissipative and purity gaps. Both these gaps must be 
finite in the bulk in order for equation (51) to hold and — as argued in section 3.2 — for the 
bulk to have a well-defined topological nature with associated topological invariants that 
are quantized. Any of these two gaps can in principle close at an edge, which gives rise 
to the variety of possibilities allowed by equation (51). In the following section, however, 
we focus on the case where the purity gap remains open at the edge, thus forbidding 
the appearance of intrinsic Majorana zero-purity modes. In that case, the topological 
properties of the steady state are equivalent to that of a pure state — as argued in 
section 3.2 — and the bulk-edge correspondence describes the behavior of a single gap 
at the edge (namely, the dissipative gap) exactly as in the Hamiltonian setting. Prom a 
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Figure 2. Left: Translation-invariant dissipative dynamics terminated at the edge 
in the case of "cross-shaped" Lindblad operators acting on a subset of 5 sites (as in 
the exphcit example of section 8). Lindblad operators acting within the system (in 
blue) are taken into account into the dissipative dynamics, whereas Lindblad operators 
requiring truncation at the edge (in red) are not considered. Right: Lindblad operator 
Li associated with a lattice site i that coincides with its center of symmetry Si, and 
example of an edge whose orientation is according to a vector e_L (see main text). 



topological point of view, all steady-state properties are the same as if the steady state 
were completely pure. We therefore restrict ourselves, without loss of generality, to pure 
(Gaussian) steady states. In addition we assume translation invariance in the following 
in order to make the explicit construction of Majorana zero-damping modes analytically 
tractable. 

We consider a (i-dimensional lattice system with a (d — l)-dimensional physical edge 
and assume that the system evolves under a translation-invariant dissipative dynamics 
consisting of a periodic repetition, on the lattice, of a single quasi-local dissipative 
process (or Lindblad operator) everywhere in the bulk and as close as possible to the 
edge (see figure 2), so that every lattice site i of the bulk becomes associated with a 
Lindblad operator Li whose form is independent of i. In order to ensure that the steady 
state is pure we further assume that the Lindblad operators form a set of anticommuting 
operators. For simplicity (and in order to incorporate the typical physical constraints 
discussed in section 4), we moreover restrict ourselves to Lindblad operators Lj that 
possess a center of symmetry, which we denote as Si (the index i used to distinguish 
Lindblad operators then corresponds, by convention, to the index of the lattice site 
located closest to Si *). We define a position vector corresponding to Si and similarly 
define the position of the lattice sites j as rj. In the translation-invariant setting defined 

*If there exists many such points, an arbitrary convention can be chosen. 
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above, the Lindblad operators then take the generic form 

Lj = ^ Uj-i aj + Vj^i a], (52) 

where denotes the subset of sites j onto which Lj acts in a non-trivial way (see 
figure 2), and uj-i, Vj-i, and aj are shorthand notations for u{rj — Rj), v{rj — Rj), 
and a{rj), respectively *. As argued in section 2.7 above, a necessary and sufficient 
condition for a Majorana mode 7 to correspond to a Majorana zero-damping mode of 
the dissipative dynamics is {Lj,7} = (for all i). Using this condition as well as the 
translation-invariant form of the Lindblad operators, one can show (see Appendix A. 3) 
that any Majorana zero-damping mode 7 must take the generic form 

d 

^^j^e'^i/^ J2 (^eir(^e,r-..(^e,ra(r, + 5]m„e„) + /i.c. , (53) 

{mi,m2, n=l 
...,m^} 

where > is a normalization factor, 0j e [0, 27r) a phase, {en}n=i,2,...,d a set 
of primitive vectors associated with the d-dimensional Bravais lattice on which the 
system is defined, and {mi, 777.2, ■ ■ ■ , m^} a set of integers defined such that the vectors 
Fj -|- Yl'^=i '^n^n span the positions of all sites in the system. Most importantly, 
{/3e„}n=i,2,...,d IS a sct of real factors determining the increase of the "spatial wave 
function" corresponding to 7 in each of the directions defined by the primitive vectors 
e„. Note that = (/3j_j)~^ is satisfied owing to the translation invariance discussed 
above {Pj-i being a shorthand notation for /3(rj — r^)) f. 

Equation (53) states that 7 — or, more precisely, its spatial "wave function" — is 
completely delocahzed when |/?e„| = 1 for all n. In the less pathological case where 
|^e„| 7^ 1 for some 77, 7 grows exponentially in the direction (7„e„ with (7„ = ±1 defined 
such that \Panen\ > 1 (or, equivalently, such that |/5_(^„e„| = |/5o-„e„r^ < !)■ In that 
case, 7 can be normalized if and only if the system is finite in the direction of "steepest" 
exponential increase, i.e., as long as the edge of the system is perpendicular to the 
direction defined by the unit vector oc Xln=i ^^S (l/^o-nen D^nGnj as shown in figure 2. 
The corresponding Majorana zero-damping mode 7 is then exponentially localized along 
that edge, decaying in every direction — (T„e„ away from the edge (into the bulk) on a 
characteristic length scale ^„ = 1/log (|/3o-„e„|)- 

In general, the existence of Majorana zero-damping modes 7 of the form (53) can 
be assessed from the knowledge of the Lindblad operators by looking for solutions of 
the equation 

0= [c-^'^*i.,_,(^,_,-^._J+^;,_,(^,_, + /3._j] (54) 

*Note that the form of the coefBcients and Vj-i ensures that has the same structure 
independently of i, as required by translational symmetry. 

fNote that is defined with respect to the site i associated with Lj, whereas Uj_j and Vj-i are 
defined with respect to the center of symmetry of ij. 
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which can be expressed solely in terms of {Pe„}n=i,2,...,d using the relations 

(i) Pk-i = Pk-jPj-i (for any triple of indices k)); 

(ii) A-i = A = i, 

where Uj-i and vj-i are fixed coefficients defining the Lindblad operators (see 
equation (52)) and denotes pairs of sites located symmetrically around the center 
of symmetry Si of the Lindblad operators. If there exists a set {f3e„}n=i,2,...,d of real 
factors satisfying equation (54) for some value of the phase (pi, the system supports 
at least a single Majorana zero-damping modes of the form (53). We remark that the 
phase (pi has a direct physical meaning in 2D systems where time-reversal symmetry is 
broken (symmetry class D), where it is simply defined by the orientation of the edge. 
We refer to Appendix A. 3 for further details as well as for an explicit derivation of the 
above results. 

5.3. Robustness of Majorana zero-damping modes and role of topology 

We now investigate the robustness of Majorana zero-damping modes against weak 
quasi-local dissipative perturbations in the general case where the steady state is not 
necessarily pure but is characterized by a finite bulk purity gap, so that its topological 
nature is well-defined. We will assume that the dissipative dynamics giving rise to such 
modes consists of n independent dissipative processes (described by Lindblad operators 
Li (i = 1,2, ...,n)) and will distinguish two types of perturbations: (i) weak quasi- 
local perturbations affecting each of these n dissipative processes and (ii) perturbations 
that introduce additional weak quasi-local dissipative processes (or Lindblad operators). 
Distinguishing spurious from genuine Majorana zero-damping modes as in section 2.7, we 
will demonstrate that topology does not guarantee the robustness of genuine Majorana 
zero-damping modes against perturbations, although it is crucial to be able to isolate 
such modes spatially. Robustness may additionally be ensured by geometrical properties 
of the system imposed through quantum reservoir engineering. 

Majorana modes are mathematically defined in the mode space A4 = consisting 
of real linear combinations of local Majorana basis operators Cj (j = 1, 2, . . . , 2A^), where 
N is the number of lattice sites in the system. In order for a particular Majorana mode 
7 = 7^ to be an exact zero-damping mode of the dissipative dynamics, one must have 
{Lj,7} = {L|^,7} = for all i (see section 2.7), which translates, in mode space, as 
two independent conditions {Li + L],'-/} = and {Lj — I/J,7} = for each i. Every 
Lindblad operator therefore imposes two independent constraints for Majorana zero- 
damping modes in as long as it is neither Hermitian nor antihermitian, which 
can safely be assumed, e.g., in the presence of disorder. Consequently, the dynamics 
generated by n < N dissipative processes necessarily gives rise to 2{N — n) exact 
Majorana zero-damping modes, independently of the form — quasi-local or not, with or 
without specific symmetries — of the corresponding Lindblad operators, as illustrated 
in figure 2. We emphasize, however, that such modes which trivially do not take 
part in the dynamics need not be spatially localized and isolated from each other (as 
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is required for them to exhibit non-Abehan exchange statistics; see section 6). As 
embedded in the dissipative bulk-edge correspondence, it is the topological nature of 
the bulk that plays a crucial role in isolating Majorana zero-damping modes away 
from each other, allowing to "fractionalize" the fermionic degrees of freedom of the 
system. Robust isolated Majorana zero-damping modes may therefore arise from the 
interplay between topology and the "incomplete" nature of the dissipative dynamics 
(i.e., n < A^), as we will demonstrate through explicit examples in section 8 below. 
Robustness in that case stems from the fact that the number n of dissipative processes 
(or Lindblad operators) taking part in the dynamics is a well-controlled quantity in 
typical physical implementations. Intuitively, this can be understood from the fact that 
driven-dissipative processes considered here require a finite (typically large) amount of 
energy to occur and therefore do not arise unless they arc deliberately introduced via 
quantum reservoir engineering (see section 2.3). In light of this physical constraint, 
dissipative perturbations of the type (i) discussed above are the most relevant ones to 
consider in assessing the robustness of Majorana zero-damping modes. We will focus on 
the latter. 

We remark that there is a priori no guarantee to find Majorana zero-damping 
modes in the case where the system is driven by n > A^ Lindblad operators — 
independently of the topological nature of the bulk — since, as discussed above, every 
generic dissipative process introduced in the dynamics leads to two additional constraints 
in mode space and further reduces the subspace available for Majorana zero-damping 
modes. Whereas spurious Majorana zero-damping modes can generally be removed by 
introducing additional dissipative processes — leaving no trace in the purity or damping 
spectra — genuine Majorana zero-damping modes, in contrast, must either survive or give 
rise to intrinsic Majorana zero-purity modes as dictated by the dissipative bulk-edge 
correspondence. Which of these two outcomes actually occurs cannot be determined 
from topological properties but rather depends on the dissipative boundary conditions, 
i.e., on the specific form of the edge dissipative dynamics resulting for the combined 
properties of bulk engineered dissipative dynamics and the physical edge delimiting it 
spatially. As we will demonstrate using explicit examples in sections 7,8 below, the 
dissipative dynamics can be constrained by e.g. geometric properties of the edge in 
such a way that the purity gap closes in a robust, controlled manner at that edge, 
thereby giving rise to an odd number of Majorana zero-damping modes in a phase 
whose bulk topology would naively imply the existence of an even number of such 
modes. Interesting phenomena with no Hamiltonian counterparts may therefore arise 
from Liouvillian dynamics as far as the edge physics is concerned, with intriguing effects 
such as the transformation of Abelian vortices into non-Abelian ones (see section 8). 

6. Non-Abelian statistics of dissipative Majorana modes 

A key property of Majorana modes in the Hamiltonian context is their non-Abelian 
statistics under spatial exchange. An important question is therefore whether or not 
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this property is also present in the case of dissipative Majorana modes (more precisely, 
of Majorana zero-damping modes). Here we provide a simple affirmative argument 
following reference [30]: Physically, this may be understood from the fact that Majorana 
modes in the context of topological insulators or superfiuids are not dynamical, but are 
purely static degrees of freedom. The actual bulk dynamics is irrelevant; it can be either 
unitary or dissipative. 

In order to verify non-Abelian exchange statistics, we first explain in more detail 
how, in a dissipative setting, the edge subspace of the density matrix — containing the 
Majorana modes — is isolated from the bulk subspace of the latter. We then consider 
the effect of adiabatic parameter changes of the Liouville operator or, more generally, 
of the operator generating the system dynamics. 



6.1. Dissipative isolation of the edge mode subspace 

As argued in sections 2.6, 2.7 and 5, an important prerequisite for the presence of 
stable dissipative Majorana modes is the existence of a zero-mode subspace as well as 
its isolation from the bulk. Here we formulate the conditions for such a situation in a 
general framework that is valid beyond the quadratic setting introduced in section 2.4. 
To this end, we introduce projectors p and q = 1 — p on the edge and bulk subspaces, 
respectively. A decoupled edge subspace then appears if the Lindblad operators ii are 
block diagonal in this projection, i.e., ii^pq — ii^qp — 0, with an edge block identical to 
zero, i.e., ii^pp — 0. The dissipative evolution then reads 

Ppp Ppq I _ ( ^ ~2Ppi^i^l,qq^^'11 

Pqp Pqq J \ ~ 2 ^ « ^l,gg^imPqP ^QQ [Piol 

(55) 

Here, the bulk dissipative evolution jCqq[pqq] has a Lindblad form and ppp is a constant 
of motion, by construction. Crucially, the density matrix elements pqp and ppq coupling 
these sectors damp out according to pqp — e~^*^<'99^»'«9*pgp(t = 0). In the presence 
of a dissipative gap, this process is exponentially fast. In steady state, the density 
matrix then takes a factorized form p — Pedge <8) Pbuik, with Pedge = Ppp and Pbuik = Pqq- 
Clearly, the absence of dynamics in the decoherence-free subspace does not allow for a 
controlled initialization of its occupation. This property is in complete analogy with a 
Hamiltonian ground-state scenario, where the zero-energy edge subspace is decoupled 
from Hamiltonian dynamics. The preparation and detection ideas developed in [68] for 
the ground states of fermionic atoms can be directly applied to the dissipative setting. 

Moving from a second- to a first-quantized representation (for a quadratic setting 
and using a Majorana basis as in section 2.5), we obtain the evolution 

dt I ^^'^ 1 = I ^ ~^P<1-^'}<1 j ^gg'^ 

y Tgp Tpp J y~-^qq^qp ~ {-^ qq y ^ qq} ~^ '^qq J 

showing explicitly that the fadeout of bulk-edge coherences is indeed controlled by the 
dissipative gap, i.e., by the slowest rate in Xqg. 
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6.2. Adiabatic parameter changes and braiding 

We consider a steady state with Majorana zero-damping modes whose corresponding 
eigenvectors jo;) span a non-local decoherence-free subspace described by the density 
matrix elements (pedge)a/3 — {(^\p\P) — Pa/s- The decoherence-free subspace has the 
property pap = 0. We now study the time evolution of the density matrix in a co- 
moving basis \a{t)) = U{t)\a{0)) that follows the decoherence-free subspace of the edge 
modes, i.e., that preserves the property Pa/s = 0. Without specifying the actual dynamics 
generating the physical evolution, the time evolution in that basis is given by 

dtp = ^ (^\a)pab{b\ + \a)pab{b\ + \a)pab{b\^ 

a,b 

= XI \^)i'^\^)P<^b{b\ + X \a)pab{b\ + X \a)pab{b\c){c\. (57) 

a,b,c a,b a,b,c 

Wc define the vector potential Aab = {a{(})\Wlf\b{0)) = |((a|&) - {a\b)), which is 
antisymmetric and Hermitian, by construction. Taking into account the normalization 
dt{b{t)\a{t)) = 0, we then obtain 

dtp = -i[A, p] + J2 \a)Pab{b\ - -t[H + A,p] + Y.{lipl\ - p)\ (58) 

a, 6 i 

where Pab = {0'{t)\dtp\b{t)) is the time evolution in the instantaneous basis. The 
Heisenberg commutator clearly reflects the appearance of a gauge structure [69-73] in the 
density matrix formalism. Crucially, this structure emerges independently of whether 
the physical dynamics encoded in pab — inserted in the last inequality — is unitary or 
dissipative. 

The transformation applied in the zero-mode subspace of either a Hamiltonian 
or a Liouvillian starting from the initial condition Pa/3(0) is given by pa^it) = 
{V{t)p{0)V{ty)af3, with time-ordered V{t) = T exp (-i rfrA(r)) where A{t)ai3 is the 
projection of the vector potential onto the dccohcrcncc-frcc subspace. The adiabatic 
change of the parameters is the key feature for such a state transformations to be 
realized while preserving the zero-mode subspace. Here one crucially requires the ratio 
between the rate of parameter changes encoded in A and the dissipative gap to be very 
small. Due to this separation of time scales — enabled by the non-evolving subspace — 
the decoherence-free subspace is never left. This phenomenon is sometimes referred to 
as the Quantum Zeno effect [57]. 

In the first-quantized formulation, the eigenvalues and eigenvectors of X now 
depend on time. The transformation into the instantaneous basis \a(t)) = f/(t)|ao), 
where |ao) is the initial reference basis, is now orthogonal. The vector potential is 
given by the Hermitian antisymmetric matrix Aab = ^((o(^)|J'(^)) — (a(^)l^(^)))) and the 
correlation matrix evolves in this basis according to {h is defined above equation (13)) 

dtT=[h + A,T]-{X,T} + Y. (59) 

Starting form this basic understanding, one can construct adiabatic local parameter 
changes of the Lindblad operators to perform elementary dissipative Majorana 
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moves [30]. Such procedure can then be apphed sequentially in order to exchange two 
particular modes, keeping them sufficiently far apart from each other during the process 
(in networks of ID wires, such exchange can be made using a T-junction geometry [16]). 
The unitary braiding matrix describing the complete process is Bij = exp (|7j7j) for two 
Majorana modes Since [Bij, Bj^] 7^ for i 7^ j, non-Abelian statistics is obtained. 



7. Illustrative examples in ID 

In this section, we apply our results to several dissipative models in one dimension. In 
particular, we investigate the interplay between criticality and purity at the onset of non- 
equilibrium topological phase transitions, exemplifying the results of section 3.3. We 
focus on ID lattice systems of spinless fermions with translational symmetry evolving 
under a dissipative dynamics described by Lindblad operators L„ = ^^^^n-mCtm + 
, where a|„ and are creation and annihilation operators associated with a 



'^n—mOj'mi vvj.h^j.<^ "'m 



lattice site m and I'n-m a-nd Un-m are translation-invariant functions. The system is then 
most conveniently described in momentum space with Fourier-transformed Lindblad 
operators Lk — Ylim ^^^^^m of the form 

Lk = Vka\ - Uktt^k- (60) 

Here we assume that 

vl = v-k and ul = U-k, (61) 

as will be satisfied in the examples examined below. The dissipative dynamics obviously 
occurs in decoupled momentum sectors corresponding to the fermionic modes a|. and a_fc. 
Defining a corresponding basis of Majorana operators (c2fc-i, C2fc, C2(_fc)-i, C2(_fc)), the 
matrices Xk and Yk describing the dynamics in each of these sectors take the following 
form (see section 2.5): 

Y _n( iW? + \vk\^)h 2Re{ulVk)(T, \ 

M 2Re{ulvk)a, {\uk\' + \vk\')h i' ^ ^ 



Y =a{ + 2lm{ulvk)(J, I 



where (T^=x,y,z denotes the usual Pauli matrices. The system exhibits critical behavior if 

at least one of the eigenvalues A[,^'^^ = ±1^^ ± Vk\^ of Xk vanishes for some /c, i.e., if the 
dissipative gap closes. Note that in that case the matrix Yk identically vanishes, since 
its eigenvalues arc given by ±(|-u^ + '^fc|)(l'"fc — Vk\)- 

All matrices Xk of the form (62) can be diagonalized through the same unitary 
transformation. This allows us to obtain the steady-state correlation matrix in a generic 
form (solving the steady-state equation {Xk^Vk} = Yk, see section 2.5); namely, 

P 1 f {\vk\^ - \uk\'^)icTy 2lm{ulVk)a^ . 

^ \vk\'^ + \uk\'^ \ -2lm{ulvk)az {\vk\'^ -\uk\'^)i(^y 
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The symmetry class to which the steady state belongs can be readily determined from 
this expression. We recall that one can associate a fictitious free-fermion Hamiltonian 
Hr = i j TijCjCj with the steady-state correlation matrix F (sec section 3.2). Here 
we write Hr = 2 J^k ^I^fe^it in a momentum-space Nambu representation with 
= {al, a-k) and obtain 

/a \vk\'-\uk\' 2ilmiulvk)\ . . 

[^l -C-k J Nl \ -2iImK^,) \uk\' - \vk? ) ' ^ ^ 

where — \vk\'^ + \uk\'^. Equation (61) implies that ^-k — and — A^, and one 
can easily verify that time-reversal and particle-hole symmetries are satisfied as T-CLf, — 
+'Hk and crx'H*_k<7x = —Tik, respectively (see equations (46) and (47)). Consequently, 
all (Gaussian) steady states resulting from a translation-invariant dissipative dynamics 
described by matrices Xk and Yk of the form (62) and (63) have chiral symmetry and 
accordingly belong to the symmetry class BDI of Altland and Zirnbauer. 

In what follows, we present three examples of dissipative systems belonging to the 
symmetry class BDI exhibiting a topological phase transition due to the change of some 
external parameter k. We first consider an example where the dissipative dynamics gives 
rise to a pure steady state in the whole parameter range of k and the system exhibits 
critical behavior at the point where the topological phase transition occurs. We then 
discuss an example in which the topological phase transition is also driven by criticality 
but accompanied by the closure of the purity gap. Finally, we provide a third example in 
which the topological phase transition docs not lead to any critical behavior. The three 
examples presented below therefore illustrate the three possibilities for non-equilibrium 
topological phase transitions identified in section 3.3. 

We illustrate our results by plotting, as a function of k, the dissipative gap 
(given by the minimum eigenvalue of X) and the purity gap which we define here as 
Ap = minfc{Tr(r^)/4}. We also consider finite systems with open boundary conditions 
but translation-invariant Lindblad operators and investigate the behavior of the edge 
zero-damping modes in the vicinity of the topological phase transition. 



7.1. Example 1: topological phase transition of a pure state with criticality 



We consider translation-invariant Lindblad operators L„ that act on three consecutive 
sites: 

1 



(66) 



where k is a real parameter controlling the coupling to the central site n *. In momentum 
space, these Lindblad operators take the form Lk — Vka\ — Uka-k with f ^ = k -|- 2 cos k 
and Uk — — 2isin/c (note that equation (61) is satisfied). The corresponding matrix Xk 
is diagonal, namely, Xk — -\- 2«;^ -|- 8kcos/c)/(4 -|- k^)1, and one can readily see that 

*Note that the overall multiplicative factor is introduced for convenience and does not affect the 
physical properties of the system. 
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Figure 3. Topological phase transition driven by criticality. (a) The system driven by 
the translation-invariant Lindblad operators of equation (66) exhibits critical behavior 
at K = ±2, as indicated by the closure of the dissipative gap Ad. At these points, the 
system undergoes a transition from a topologically trivial {viu = 0) to a topologically 
ordered (i^id = 2) phase. Since the purity gap Ap = 1, the steady state of 
the dissipative dynamics is always pure, (b) Divergence of the localization length 
associated with one of the two Majorana zero-damping modes found at the edges upon 
approaching the topological phase transition at k — 2. 



the dissipative gap closes at k = ±2, giving rise to critical behavior around these points. 
The steady-state correlation matrix can be written explicitly as 

^ / (k^ + cos k + Acqs 2k)iay 4:{smk + Ksin2k)az 

^ \ — 4(sin A; -|- Ksin2A;)cr2 (k^ -|- 4k cos A; -|- 4 cos 2A;)i(T, 

with Qk = 1 / {4 + K,^ + cos k) and one can easily verify that = —I for all k, meaning 
that the steady state is always pure. 

Expressing the momentum-space correlation matrix in the form = i(nfc ■ cr) with 
rifc G similarly as in section 3, we obtain = gk{0, 4K(sin k + sin 2k), + 4k cos k + 
4cos2/c). As expected, this vector is non-zero for all k in the whole parameter range 
of K where the dissipative gap is finite, i.e., for |k| ^2. In that case one can "spectrally 
flatten" the steady-state correlation matrix by normalizing the vector to unity 
in order to calculate the corresponding winding number topological invariant vid (see 
equation (35) and discussion above). Here we find z/id = 2 for |k| < 2 and z/id = for 
\k\ > 2 (for an explicit calculation, see our previous work [31]). We therefore identify a 
topological phase transition at \k\ =2 where the dissipative gap closes and the system 
becomes critical. These results are illustrated in figure 3(a). 

We now examine the case of a finite system (or "chain") of N sites with open 
boundary conditions. We first notice that only N — 2 three-site Lindblad operators 
of the form (66) can be "placed" on such an open chain. Hence, according to our 
discussion of section 5.3, four exact Majorana zero-damping modes (two at each edge, 
by symmetry) must be present independently of the topological properties of the system 
(in particular, for any value of the parameter k). Obviously, this number which remains 
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constant across the topological phase transition cannot be used to identify the latter. 
However, the topological phase transition can be revealed in two other ways: (i) 
One can remove spurious Majorana zero-damping modes with no topological origin 
by introducing additional arbitrary Lindblad operators acting at the edges; assuming 
that these extra Lindblad operators preserve chiral symmetry *, the system must then 
exhibit a total of i^iD = 2 Majorana zero-damping and/or zero-purity modes in the 
parameter range \k\ < 2, as dictated by the dissipative bulk-edge correspondence (see 
section 5.1). (ii) One can examine the behavior of the localization length associated with 
the Majorana zero-damping modes as the parameter k is tuned across the topological 
phase transition; owing to criticality, the characteristic length scale associated with at 
least one of the two modes found at an edge must diverge at the phase transition point 
\k\ = 2. 

Here we follow strategy (ii): Since the Lindblad operators L„ are translation- 
invariant and have a symmetric form around their central site n f, we use the results of 
section 5.2 and assess the existence of Majorana zero-damping modes that decay into the 
bulk by solving equation (54). Choosing 0j = 0, the latter equation takes the explicit 
form — K + 2/3ej, (e^; being a unit vector joining neighboring sites of the chain) and we 
readily obtain (3^^ — —k/2, showing that a least one of the two Majorana zero-damping 
modes at the edge must decay (exponentially) into the bulk when < 2 (as required 
by \PeJ < 1) on a characteristic length scale ^ = —1/ \og{\l3eJ) = — 1/ log (|/t|/2). 
As expected, ^ diverges at \k\ = 2 and thus reveals the topological phase transition. 
Figure 3(b) illustrates numerical results corroborating this behavior. 

7.2. Examples 2 and 3: topological phase transition of a mixed state with and without 
criticality 

We now consider a dissipative dynamics that involves two types of translation-invariant 
Lindblad operators: 



acting on neighboring and next-to-nearest neighboring sites, respectively. The Lindblad 
operators Ln^ generate a dissipative dynamics that drives the system into a pure steady 
state corresponding to the ground state of a topologically non-trivial Kitaev chain [15], as 
demonstrated in our previous work [30]; two Majorana zero-damping modes are found in 
that case (one at each edge) . The Lindblad operators Ln^ describe two decoupled Kitaev 
chains, with four Majorana zero-damping modes (two at each edge). Below we discuss 
two scenarios: (i) We first assume that both dissipative processes L^'' and L^^ occur 

*Note that time-reversal symmetry only must be ensured since particle-hole symmetry is 
automatically satisfied in our dissipative framework. 

f Namely, the creation and annihilations parts of -L„ (see equation (66)) have s-wave and p-wave 
symmetry, respectively. 




(68) 



(69) 




Figure 4. Generalizations of the dissipative Kitaev chain. By positioning the sites 
of the ID chain (grey disks) in a zigzag geometry, three different types of dissipative 
processes can be engineered (see section 2.3) from auxihary sites (blue circles): (a) 
Lindblad operators Ln"^ driving the system into the ground state of the Kitaev chain. 

(2) 

(b) Lindblad operators L„ driving the system into two decoupled Kitaev chains, (c) 

By moving the auxiliary sites relative to the physical ones in the zigzag geometry, one 

can engineer a coherent superposition L„ cx lI^'' +KLn^ of the Lindblad operators L^n"^ 

and Ln . (d) Furthermore, by doubling the number of auxiliary sites (introducing the 

yellow ones), one can engineer a dissipative dynamics £ oc £i + corresponding to 

two "competing" Liouvillians Ci and £2 associated with Lindblad operators L^'n^ and 
(2) 

Ln , respectively. 



coherently, so that the relevant Lindblad operators are L„ = {Ln^ + kLu'^) / {2{k'^ + k + 1)) 
with K G M. (ii) We then consider the case where both dissipative processes "compete" 
against each other, so that the relevant Liouvillian to describe the dissipative dynamics 
is £ = (£1 + /t£2)/(l + 1^), where Co, denotes the Liouvillian corresponding to a single 
Lindblad operator L^n^ and k > 0. In practice, Ln"* and Ln'' as well as their combinations 
(i) and (ii) can be realized in a "zigzag" geometry, e.g., as depicted in Fig. 4. 



In momentum space, the Lindblad operators take the generic form = v^^^a]. 



i"^a_fc (a = 1,2) with Vk = cos{ka/2) and Uk = sm{ka/2) (note that 

equation (61) is again satisfied). We first examine case (i), in which the eigenvalues of 
the matrix take the explicit form 

(l + /t)2 , l + 2Kcosk + K,^ , , 

showing that the dissipative gap closes at = ±1, and the steady-state correlation 
matrix reads 

I (cos A; + K cos 2A;)icrj; — (sin A; + /tsin2fc)cr2 , 

(sin /c + Ksin2A;)o"2 (cos A; + k cos 2/c)icrj^ 
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Figure 5. Topological phase transition of a mixed state with criticality. (a) The 
system driven by translation-invariant Lindblad operators of the form L„ cx L^n^ +KL^n'^ 
exhibits critical behavior at /c = ±1, as indicated by the closure of the dissipative gap 
Ad- At these points, the system undergoes a topological phase transition from a 
state characterized by v^d = 1 to a state with i^iD = 2. The steady state is mixed 
(unless K = or K — > zbcx)) and the phase transition is crucially accompanied by the 
closure of the purity gap Ap. (b) Divergence of the localization length associated with 
the Majorana zero-damping mode found at each of the edges upon approaching the 
topological phase transition at k — 1. 



with Qk = {1 + k) / {1 + K + + KCOS k) and eigenvalues ±(1 + + + 2k cos k/{l + 
K + + Kcosfc), indicating that the steady state is pure for k = and k — )■ ±oo 
and mixed otherwise (as expected from the fact that k = (k — )■ ±oo) corresponds to 
the case where Ln'' (respectively L^n^) alone generates the dynamics). The vector 
defined so that = i(nfc ■ cr) is non-zero for all k provided that k ^ ±1 and takes 
the explicit form = (7^(0, — (sin A; -|- /€sin2/c), — (cosfc -|- k cos 2k)). The corresponding 
winding number topological invariant is z/id = 1 for < 1 and Uij) = 2 for |k| > 1, in 
accordance with the fact that = 1 for a single topologically non-trivial Kitaev chain. 
These results, which are depicted in figure 5(a), illustrate the possibility of a topological 
phase transition driven by criticality and by the closure of the purity gap. 

To conclude our investigation of case (i), we consider a finite chain of sites 
with open boundary conditions. Since only N — 2 Lindblad operators of the form 
can be applied on such a chain, we again obtain four exact Majorana 
zero-damping modes that are either genuine or spurious. A numerical analysis below 
the topological phase transition point k, = 1 reveals that two of these modes (one at each 
edge) decay exponentially into the bulk on a characteristic length scale ^ = — l/log(K) 
which diverges at expected (see figure 5(b)). 

We now turn to the investigation of case (ii) in which the two dissipative processes 
Ln^ and Ln^ compete against each other. In that scenario, one can readily verify that 
Xfc = 21, such that the damping spectrum is gapped (and "flat") for all k. The steady- 
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Figure 6. Topological phase transition of a mixed state without criticality. The 

system evolving under a dissipative dynamics £ oc £1 + generated by two 

"competing" Liouvillians L\ and £2 corresponding to the Lindblad operators L^n^ 
(2) 

and Ln , respectively, exhibits a topological phase transition at k = 1 from a state 
with winding number i^id = 1 to a state with winding number v\y) — 2 that results 
from the closure of the purity gap Ap only. The system never becomes critical, since 
the dissipative gap Ad remains constant in the whole parameter range of n. 



state correlation matrix is then given by 

^ ( {cosk + Kcos2k)iay — (sin/c + /tsin2A;)cr2 
(sin /c + K sin 2A;)cr^ (cosA; + Kcos2A;)icry 

with g = 1/(1 + k) and eigenvalues ±Vl + + 2fi;cos A;/(l + k), implying that the 
steady state is pure if and only if k = or /t — i- 00. The corresponding vector is 
non-zero (for all k) except at k = 1 in which case Tk=±TT = 0. For k 7^ 1, one finds 
rifc = g{0, — (sinfc + /tsin2fc), — (cos + k cos2fc)) and the corresponding winding number 
topological invariant is z/id = 1 for k < 1 and uid = 2 for k > 1. These results, which 
are presented in figure 6, demonstrate that case (ii) provides an example of a topological 
phase transition that is not driven by criticality but rather by the closure of the purity 
gap only. 

Examining case (ii) on a finite chain with open boundary conditions, we find four 
Majorana zero-damping modes (two at each edge) in the whole parameter range of k. In 
accordance with the dissipative bulk-edge correspondence (see section 5.1), all of these 
modes must be genuine for k, > 1 where i^id = 2, while two of them (one at each edge) 
must be spurious for k < 1 where z/id = 1. Crucially, the spatial wave function of these 
four modes does not exhibit any critical behavior * upon approaching the topological 
phase transition point k, = 1. 




*In fact, it does not change at all in this particular example. 
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8. Illustrative example in 2D 

In the previous section, we have exemphfied the phenomenology of Majorana zero- 
damping modes at physical edges in the context of ID systems. Here we focus instead on 
the physics that arises when topological defects are introduced in the bulk of the system, 
away from physical edges. We consider topological defects in the form of dissipative 
vortices in infinite 2D systems in order to illustrate such physics in the simplest possible 
scenario. 

8.1. Dissipative vortices 

Dissipative vortices are most conveniently defined starting from an infinite 2D lattice 
system evolving under a dissipative dynamics generated by translation-invariant 
Lindblad operators Lj = J2j '^j-i aj + Vj^i (using notations as in section 5.2 above; see 
equation (52), in particular). We assume that the corresponding Liouvillian exhibits a 
finite dissipative gap and that its steady state is characterized by a finite purity gap, so 
that the topological nature of the system is well-defined (i.e., the steady state belongs to 
a specific topological class and is characterized by quantized topological invariants); in 
particular, the dissipative bulk-edge correspondence introduced in section 5 is satisfied *. 
We then introduce a dissipative vortex at a position defined by the vector Rq f by 
modifying the translation-invariant coefficients = Uij defining the annihilation part 
of the Lindblad operators in the following way: 

Uij^Uijf{rj)e-''^^, (73) 

where (r^, ipj) are polar coordinates defining the position of each site j with respect to 
Ro, /(r) is a real and positive function that vanishes as r — > and reaches a constant 
value as r ^ S {S being a characteristic length scale associated with the dissipative 
vortex, defining the vortex core), and e~^'^^ describes the vortex phase winding £ times 
around the origin Rq {£ defining the so-called vorticity). While vortices exhibit similar 
properties in Hamiltonian systems, the specific form of a dissipative vortex defined 
above is motivated by its natural realization in typical implementation schemes with 
cold atoms based on optical vortex imprinting (we refer to our previous work [31] for 
further details). 

We remark that the qualitative features of a dissipative vortex do not depend on 
the explicit form of the function /(r) describing its core. Of crucial importance is the 
fact that /(r) vanishes as r — > 0, so that the vortex core can be identified with a region 
of space in which the system is topologically equivalent to the vacuum, with all sites 
essentially occupied J. If the topological nature of the bulk is non-trivial, a small circular 

*One can assume, without loss of generality, that the Lindblad operators form a complete set of 
anticommuting operators, so that the system is driven into a pure steady state independently of the 
initial conditions (see discussion of section 3.1). 

fNote that this position need not coincide with a lattice site. 

|This can be understood from the fact that the Lindblad operators acting in the vortex core have 
an annihilation part that essentially vanishes, namely, Lj = Cj + with A, « 0. 
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edge (or domain wall) of radius ~ S must then appear around the vortex core, separating 
this vacuum-like region from the surrounding bulk which is potentially topologically non- 
trivial. Remembering the dissipative bulk-edge correspondence discussed in section 5, 
one then naively expects Majorana zero-damping modes and/or intrinsic Majorana zero- 
purity modes to appear at this domain wall. We argue below that the possibility of 
having such modes, however, additionally and crucially depends on the phase winding 
i of the dissipative vortex, as expected by analogy to the Hamiltonian setting. 

The crucial role of the vortex phase winding £ in prohibiting or allowing the 
existence of Majorana zero modes in the vortex core can be understood in analogy 
to the Hamiltonian setting. In the Hamiltonian context, circular domain walls typically 
trap low-energy modes with quantized angular momentum m and a corresponding 
energy E m [12]. When a flux of the f/(l) gauge field associated with the fcrmions 
corresponding to n quanta of angular momentum is inserted in the region enclosed by 
the domain wall, the angular momentum of these low-energy modes shifts to values m+n 
and the energy of the latter shifts accordingly, thereby prohibiting or allowing for zero- 
energy modes. We argue that the same argument applies here in the case of a dissipative 
vortex, namely: the phase winding £ defines a fiux tt^ threading the vortex core and shifts 
the eigenvalues corresponding to low- damping or low-purity modes bound to the vortex 
core (in the damping or purity spectrum, respectively). This can be understood from the 
fact that the Hamiltonian phenomenology discussed above automatically applies to the 
parent Hamiltonian associated with the dissipative dynamics, whose zero modes either 
correspond to Majorana zero-damping modes or to intrinsic Majorana zero-purity modes 
(see discussion of section 2.7). In analogy with the Hamiltonian setting (see, e.g., [12]), 
we then expect Majorana zero-damping modes and/or intrinsic Majorana zero-purity 
modes to appear in dissipative vortices with odd vorticity £ only, as we corroborate 
below in an explicit model. 

We remark that the physics associated with a dissipative vortex remains 
qualitatively the same upon insertion or removal of 2ti fiuxes in the vortex core since 
such modifications can be seen as gauge transformations [12]. The physical properties 
of a dissipative vortex are therefore determined by the vorticity I modulo 2, and one can 
restrict oneself to £ = or £ = 1, without loss of generality. Since the case of a dissipative 
vortex with £ = is physically equivalent to that of a small circular physical edge, we 
will focus exclusively on dissipative vortices with £—1. Remarkably, the physics arising 
in a single dissipative vortex with a tt flux {(■ — 1) in an infinite system on the plane 
can be shown to be equivalent that occurring at the edge of a semi-infinite system with 
a cylinder geometry and no fiux (see figure 7 and our previous work [31] for further 
details). This equivalence — which will be useful in the sections below — can intuitively 
be understood from the fact that the tt fiux of the U{1) gauge field introduced by a 
vortex with £ = 1 naturally arises in cylinder geometry due to the (extrinsic) curvature 
of the latter. 
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Figure 7. Left: Semi-infinite plane "wrapped up" into a semi-infinite cylinder in the 
Gy direction (shown in red). Right: Topological equivalence between the semi-infinite 
cylinder geometry and an infinite planar geometry with a single dissipative vortex 
"puncturing" the plane. 



8.2. Explicit 2D model and hulk properties 

We now consider, in the framework defined above, an explicit model defined on a square 
lattice with primitive vectors and e^, with translation- invariant Lindblad operators 
Li = X]j "^j-i + "^j-i ^] whose only non-zero coefficients are 

M±e, = ±1; M±e, = ±i, (74) 

where /3 is a real parameter which can be used to tune the system across phase 
transitions. Every site i therefore has an associated Lindblad operator Lj of the form 

Li = Cj + A = /3a\ + {al^^ + + al^^ + alj 

+ (aj+e, + icii+ey - Cli-e^ " ICtj-eJ, (75) 

where the creation part Cj is isotropic or s-wave symmetric with respect to the 
central site i (as naturally arises in typical implementation schemes; see section 2.3) 
and the annihilation part Ai has a p-wave symmetry, i.e., Ai = (Vx + iVy)aj with 
Vxtti = flj+e^ ~ o,i~ex- momcutum space, these translation-invariant Lindblad 
operators take the form = u^a^ + Vka^_^, where 

Uk = 2i(sin (kx) + i sin (fcy)); 

= (3 + 2(cos (kx) + cos (ky)). (76) 

Since the above functions satisfy Wk^k = — "U.kt'-k; the Lindblad operators form 
a complete set of anticommuting operators and the system is driven into a pure 
Gaussian steady state corresponding to the ground state of the parent Hamiltonian 
-f^parcnt = J2k ^k^^ (^^^ equatiou (39) and discussion thereof). The damping spectrum, 
which coincides with the spectrum of -^parent in that case, is then given by the norm 

ATk = y^{4,Lk} = Vl^kP + KP. (77) 
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Remembering the explicit form of the functions and f k defined in equation (76) above, 
one can easily verify that the corresponding dissipativc gap closes at the parameter 
values f3 = 0, ±4. We thus expect the existence of four distinct — possibly topological — 
phases depending on the value of /3. 

We remark that the Lindblad operators are not invariant under time-reversal 
symmetry due to the p-wave symmetry embedded in their annihilation part. The system 
thus belongs to the symmetry class D of Altland and Zirnbauer (see section 3.2) and 
is characterized by the Chern number topological invariant defined in equation (32). 
One can easily verify that, in accordance with the discussion of section 4.2, the Chern 
number vanishes in the whole ^ parameter range where the system exhibits a finite 
dissipative gap (we refer to our previous work [31] for an explicit proof). As dictated by 
the dissipative bulk-edge correspondence introduced in section 5, we therefore naively do 
not expect to find Majorana zero-damping modes in the system if edges are introduced. 
Despite this conclusion, let us now try to construct such modes explicitly in the 
translation-invariant setting of section 5.2 anyway. 

8.3. Construction of Majorana zero-damping modes 

The possibility of having Majorana zero-damping modes at a physical edge when the 
dissipative dynamics is generated by translation-invariant Lindblad operators as in 
section 5.2 above can be assessed from the exphcit form of the Lindblad operators 
(see equation (74)) by looking for a solution of equation (54). Choosing 0j = (which 
corresponds to choosing a particular orientation of the edge; see section 5.2), the real 
and imaginary parts of equation (54) here take the explicit form 

= ^ + 2/3e.+/3e, + (^eJ-\ 
= ^e,-(^ej-^ 

One thus finds (3ey = ±1 and f3e^ — — (^/2 ± 1), which implies the existence of a 
solution with l^g^l < 1 if and only if < \/3\ < 4. Consequently, there exists at least one 
Majorana zero-damping mode 7 of the form defined by equation (53) when < \/3\ < 4 if 
the system possesses an edge in the — e^; direction *, in which case 7 is uniformly spread 
along the edge (since \Pey\ = 1) and exponentially localized in the -l-e^; direction away 
from the edge on a characteristic length scale ^ = — 1/ log (|/9e^|) = — 1/ log (1/3/2 ± 1|). 
The semi-infinite planar geometry of the system in that case is depicted in figure 7(a). 
Remarkably, the special points /3 = 0, ±4 where the dissipative gap closes coincide 
with those at which \/3e^\ = 1, i.e., at which the localization length associated with 7 
diverges. This suggests the possible existence of topological phase transitions occurring 
at the gap-closing points =0, ±4 which are not identified by the (vanishing) Chern 
number topological invariant. 

*Note that a different orientation of the edge would be obtained for <pi ^ 0. 
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8.4- Topological origin 

We have demonstrated above that the system can support Majorana zero-damping 
modes at an edge (in the parameter range < \/3\ < 4) despite its vanishing 
Chern number. We argue below that this apparent contradiction stems from the fact 
that we have assumed translational symmetry, and estabhsh the topological origin of 
Majorana zero-damping modes in the above translation-invariant setting. It is clear 
that translational symmetry is not robust against disorder. However, assuming such a 
symmetry will allow us to demonstrate, in the simple model above, a key mechanism 
with no Hamiltonian counterpart through which spatially isolated Majorana zero- 
damping modes can be obtained in a topological phase characterized by an even integer 
topological invariant, in contradiction with the Hamiltonian bulk-edge correspondence. 

We first remark that the construction of section 8.3 above remains unchanged 
if we assume that the system is finite in the ey direction, "wrapped up" into a 
cylinder (as depicted in figure 7(b)). Since such a semi- infinite cylinder geometry is 
topologically equivalent to that of an infinite plane "punctured" by a single dissipative 
vortex (figure 7(c)), Majorana zero-damping modes of the form found in section 7 must 
therefore similarly appear in the core of a dissipative vortex on the plane *. In these two 
equivalent geometries, the topological origin of Majorana zero-damping modes can be 
identified by moving to momentum space in the ey direction associated with translational 
(or rotational, in the vortex case) symmetry (see figures 7(b) and (c)). The system 
then reduces to a "stack" of semi-infinite ID wires (in the direction) associated 
with different momenta ky. In particular, the Fourier-transformed counterpart of the 
Lindblad operators defined in equation (75) take the form 

Li = KaJ + (aj+^^ + aJ_^J + (a^+e. - Oi-eJ (78) 

in the momentum sectors corresponding to ky — or tt, where i indexes the lattice 
sites in the remaining spatial direction and k — p + 2 ion: ky — and /? — 2 for 
ky — TT, respectively. In these two sectors, the system thus reduces to the ID wire 
with chiral symmetry investigated in section 7 above (up to an overall dissipation rate 
\/4+~/? which does not affect any of the topological properties of the system). Owing 
to chiral symmetry, such a system belongs to the symmetry class BDI of Altland and 
Zirnbauer and is characterized by the winding number topological invariant z/id defined 
by equation (35). As discussed in section 7 and depicted in figure 3 thereof, one finds 
i^iD = 2 for |«;| < 2 and z/id = for |«;| > 2. Discontinuities in the winding number 
at K = ±2 (or, equivalently, at ^ = 0, ±4) clearly establish the occurrence of non- 
equilibrium topological phase transitions at the dissipation gap-closing points identified 
in section 8.3. 

In summary, we have identified a non-trivial topological invariant z/id = 2 in the 
parameter range < \f3\ < 4 f by taking advantage of translation invariance to reduce 

*Note that the translation symmetry along the circumference of the cylinder translates as a 

rotational symmetry around the vortex core. 

fNote that /3 > (/3 < 0) corresponds to the ID wire in the momentum sector ky = {ky = tt). 
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the original 2D problem in semi-infinite cylinder geometry or, equivalently, in planar 
geometry with a single £ = 1 dissipative vortex, to a ID wire defined along the axis of 
the cylinder or in the radial direction away from the vortex core (e^^ direction shown 
in figures 7(b) and (c)). This provides a topological explanation as to why we were 
able to demonstrate the existence of at least one Majorana zero-damping mode for 
< 1^1 < 4 in section 8.3. In fact, invoking the dissipative bulk-edge correspondence 
of section 5.1, we can now argue that a total of i^iD = 2 Majorana zero-damping modes 
and/or Majorana zero- purity modes must be present at the edge of the semi- infinite 
system in cylinder geometry or, equivalently, in the core of the £ — 1 dissipative vortex 
on the (infinite) plane, as long as the symmetries of the system are preserved *. 

The precise number {mdamping) edge < 2 (see equation (51) and discussion thereof) 
of genuine Majorana zero-damping modes that is found at the edge of the cylinder or 
in the vortex core depends on the specific dissipative boundary conditions that appear 
there. Crucially, physically distinct boundary conditions naturally emerge in these two 
systems, although their geometries are topologically equivalent. In the case of the 
dissipative vortex, in particular, one can demonstrate explicitly that a single Majorana 
zero-damping mode is found in the vortex core (alongside with a single Majorana 
zero-purity mode such that the dissipative bulk-edge correspondence is satisfied; see 
equation (51)). The key mechanism through which the dissipative vortex isolates exactly 
one Majorana zero-damping mode despite the fact that the bulk is characterized by an 
even (integer) topological invariant (z^id) originates from the geometry of the vortex 
core alone and is therefore generic. We refer to our previous work [31] for an explicit 
proof of this statement. 

8.5. Physical properties of dissipative vortices carrying single Majorana zero-damping 
modes 

In the Hamiltonian context, tt fiux vortices carrying single Majorana zero-energy modes 
have been demonstrated to exhibit non-Abelian exchange statistics [8] . The underlying 
proof exclusively relies on the algebraic properties of Majorana modes and, crucially, 
does not depend on the dynamics giving rise to them (see section 6.2). As a result, the 
same conclusions apply in the dissipative setting of interest in this work. In particular, 
TT fiux (i.e., £ = 1) dissipative vortices carrying single Majorana zero-damping modes 
must exhibit non-Abelian exchange statistics (when moved around through adiabatic 
parameter changes as discussed in section 6.2). In light of the mechanism illustrated 
in the above explicit model, we thus conclude that vortices can have, in our dissipative 
setting, non-Abelian exchange statistics in a phase characterized by an even integer 
topological invariant, in stark contrast to what can occur in the Hamiltonian context. 

*More precisely, translational symmetry along the circumference of the cylinder or, equivalently, 
rotational symmetry around the vortex core must be preserved, as well as the chiral symmetry of the 
dimensionally reduced problem (ID wire). In the presence of disorder, one expects Majorana zero- 
damping modes to become "softer" , acquiring a small but finite damping rate. 
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Figure 8. Numerical results for a configuration of two i — 1 dissipative vortices on 
a square lattice of 35 x 35 sites with unit spacing, (a) Generic form of the Majorana 
zero-damping modes localized in each of the vortex cores (here for /3 — 2 and a vortex- 
vortex separation d = 16), and (b) exponential behavior of their residual damping 
rate as a function of the vortex-vortex separation d for /3 = 2.8 to 3.2 by steps of 
0.1 (from bottom to top). Note that, due to the finite size of the system, vortices 
"interact" more strongly (via the edge) as they come closer to the edge, giving rise 
to the symmetric behavior that can be seen around d ~ 16. (c) and (d); Low-lying 
part of the damping (respectively purity) spectrum for (3 — 2 featuring a gap with two 
quasi-zero (^ 10^^^) eigenvalues (indicated by a red arrow) corresponding to Majorana 
zero-damping (zero-purity) modes trapped in each of the vortex cores. All results were 
obtained for vanishingly small vortex cores, i.e., with /(r) = 1 everywhere except at 
r = where it vanishes (see equation (73)). 

In addition to their non-Abelian statistics, vortices carrying single Majorana zero 
modes behave, in our dissipative setting, in a fully analogous way as in the Hamiltonian 
context. In particular, two such vortices "interact" as the distance d between the latter 
becomes small, acquiring an exponentially small residual damping rate ~ e~'^/^, where 
^ is the localization length associated with the Majorana zero-damping mode that they 
carry. We have verified such phenomenology through numerical simulations. Figure 8 
summarizes the results obtained in the case of a finite system on the plane with two 
"interacting" £ = 1 dissipative vortices. 
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9. Conclusions 

In this paper we have provided a discussion of non-equihbrium superfluid topological 
states generated by means of engineered dissipation, highlighting analogies and 
differences to the more conventional Hamiltonian setting. Similarities mainly stem 
from the fact that topological properties are properties of the state of the system — 
encoded in the static correlations — ^with some care to be taken regarding the purity of the 
state. In particular, identifying the correlation matrix as a tool for the symmetry-based 
topological classification of arbitrary Gaussian states analogous to a first- quantized 
quadratic Hamiltonian, it is clear that the classification of topological states according 
to the symmetry classes of Altland and Zirnbauer extends to a general non-equilibrium 
steady-state situation, highlighting the universality of this classification. Differences, on 
the other hand, mostly stem from the fact that the spectral or dynamical properties 
of the system and the properties of the steady state are not as tightly related as 
in the Hamiltonian ground-state scenario. The interplay between the well-known 
characteristics of topological states and the new elements brought in by the dissipative 
nature of the dynamics can give rise to effects with no Hamiltonian counterpart. 
Examples are topological phase transitions which do not require the dissipative spectral 
gap to close, or the trapping of unpaired Majorana modes in a bulk with vanishing 
Chern number. 

We have presented in this work a detailed discussion of dissipatively induced 
topological superfluid states in ID and 2D systems which is a flrst step towards a 
complete picture of non-equilibrium topological order. In this respect, various directions 
remain to be explored in future work: 

Prom the viewpoint of dissipation engineering, we may ask ourselves whether other 
symmetry classes than the ones identified in this work can be achieved in realistic setups. 
So far we have concentrated on topological superfluids of spinless fermions. Adding spin 
degrees of freedom opens up the perspective of reaching topological insulating phases 
with potentially more robust edge modes. It will also be interesting to investigate 
whether going beyond quasi-local Lindblad operators is possible in realistic cold atom 
experiments, allowing to reach phases characterized by a non- vanishing Chern number. 

Prom a many-body perspective, an interesting direction will be to explore the nature 
of topological phases and phase transitions in more detail. This concerns, in particular, 
the interplay of dissipative and Hamiltonian dynamics, which we did not examine in 
this work. While this intentional omission is not a severe limitation in the context 
of cold atomic gases where the Hamiltonian dynamics can be suppressed, including 
the effects of such dynamics would give a more complete picture of the robustness of 
topological phases under general combined unitary and dissipative dynamics. From the 
viewpoint of critical phenomena, topological phase transitions accompanied by a closure 
of the dissipative gap offer intriguing examples of criticality in fermionic systems. Going 
beyond a mean-field approach by including the effects of Hamiltonian and dissipative 
interactions will be necessary to reach a detailed understanding of this issue. 
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Prom a quantum information point of view, it will be important to investigate 
in more detail the analogies and differences between dissipatively induced non-local 
decolierence-free subspaces and more conventional Hamiltonian edge mode subspaces. 
To answer these questions, a more systematic classification of harmful dissipative 
perturbations as well as a more detailed understanding of the dynamics of excitations on 
top of the dissipatively stabilized state will be necessary. The propagation of excitations 
or defects on top of a state which is entirely induced by dissipation can be expected to 
be vastly different from a Hamiltonian ground state. Also, manipulation and readout 
tools tailored to the dissipative setting will have to be developed. 
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Appendix A. Dissipative mean-field theory 

We show in more detail how the microscopically number conserving implementation of 
the Lindblad operators, leading to an interacting (quartic in the fermion operators) 
Liouvillian, are related to the quadratic setting analyzed in view of its topological 
properties. To this end, we first present a general relation of the number conserving 
operators to their fixed phase counterpart on the level of the dark state wave functions. 
We then show within a properly devised mean-field theory how these two settings are 
connected dynamically, with main result stating that the quadratic dynamics emerges 
naturally in the long time limit of the microscopically number conserving quartic 
dissipative evolution. 

Appendix A.l. General relation of fixed number and fixed phase Lindblad operators 

We show that for a given real space quadratic master equation with pure dark state, 
we can directly construct a number conserving quartic version and vice versa, and 
specify the respective fixed phase or fixed number exact dark state wave functions. In 
particular, in the thermodynamic limit, both descriptions become equivalent as made 
explicit below. 

We study a general Liouville operator generated by Lindblad quantum jump 
operators Jj which exhibits a pure dark state: 

C[p]= KY,{J^PJl - \{JU^.P}). Jr\D)=Q\ll, PD=\D){D\. (A.l) 
i 

We will show the following: For each fixed phase setting, specified by Lindblad operators 
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and dark state 

J. = Li = A + oiCj, \D) = |BCS, e) (X exp(Q;G'^)|vac), (A.2) 
there is a fixed number setting with * 

Ji = ii = CjA, \D) = |BCS, N) oc G^^lvac). (A.3) 

The reverse statement holds as weU. Here, as in the main text the creation (annihilation) 
parts Cj{Ai) are linear in the spinless fermion operators and obey the following 
requirements: 

(i) Translation invariance - C] = Vi^jO^^ and Ai = J2j Ui-jdj with translation 
invariant complex quasi-local position space functions Vi-j, Ui-j. The Fourier transforms 
for the creation and annihilation part are local in momentum space, 

i i 

The Fourier transformed fixed phase and fixed number Lindblad operators then read, 
respectively 

Lk^Ak + aCj, 4 = $^Ci_k^q. (A.5) 

q 

(ii) Antisymmerty - The functions have the property U-k — ii^k, f-k = l^v^, such that 
the related "wave function" (p is antisymmetric, 

(/7k = Vk/ui^ = -</7-k- (A.6) 
This implies that the set Lk form a full Dirac algebra up to a normalization, 

{Lk, Lk'} = {4, 4,} = 0, {Lk, 4} = (|«kr + \avu\^)6{k - k') (A.7) 
for all k, k'. We may introduce normalized operators via 

Zk = UkUk + avka^_^, Uk = u^/ = v^/ 

Nk = \uk\^ + \avk\^, (A.8) 

such that l^kp + \avk\'^ = 1, {Lk, L^,} = S{k — k'). 

The generator of the state is bilinear and reads in terms of the momentum space 
wave function 

Gt= J2vW-A- (A.9) 

k 

The fixed number operators satisfy [ii, N] — and thus conserve total particle 
number N — X^i^l^j- The corresponding dark states |BCS,iV) contain 2N fermions. 
Instead, in the fixed phase setting, [Lj, iV] 7^ 0. a = is an arbitrary complex 
number. Within the mean field described below, the modulus r can be related to the 

^Working with the functions |BCS,A^) for finite A'' may require an infrared momentum cutoff 
9l ^ so that N = nL'^ {n the particle filling, L the number of sites in each direction, d the spatial 
dimension), which is consistent with the thermodynamic limit N^oo,L—>oo,n = N/L'^ — )■ const. 
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average particle number. 9 has an interpretation in terms of a fixed, superfluid phase. 
This is made explicit from the expansion of the coherent state wave function 

|BCS,^) =A^exp(Q;G'^)|vac) = AT + Q;(/?kaLk4)|vac), 

AA=^',(l + |a^k^)-^/^ (A.IO) 

where A/" is a global normalization. is restricted to k such that /ca > in all primitive 
lattice directions A. 

With these preparations, wc can prove the above statement. We work in momentum 
space and start with the number conserving setting. The Lindblad operators £k = 
^qCq_j^Aq are normal ordered, such that £k|vac) — 0. The dark state property 
£i|BCS, N) = £k|BCS, N) — for all i or k is therefore equivalent to following relation 
for all k, 

[4,G'^] = X] ^^q-kMq<^qaq_k«-q = ' ^ VqU^_kip^_l,al^_j^al^ = 0. (A.ll) 

q q 

This is true if and only if ^'^^'^-^ = ^a- holds for all k, i.e., for the wave function 
(fq — Vq/uq up to a momeutum independent complex number. The latter can be 
absorbed into the constant a. 

It is easy to show the uniqueness of the dark state in the fixed phase ensemble. The 
coherent state wave function (A.IO) can be written as |BCS,^^) oc Hk -^k-^-k|vac), so 
that it is indeed the unique wave function annihilated by the full set of Dirac operators 
Lk up to normalization. Wc arc not aware of an analogous argument for the fixed 
number case (for each given particle number 2N), since the bilincars do not obey 
a simple algebra. Numerical simulations in the context of spinful Lindblad operators 
[52, 53] however suggest a unique steady dark state. 

The two wave functions |BCS, A^"), |BCS,^) are equivalent in the thermodynamic 
limit A" — > oo: We consider the relative number fluctuations of the fixed phase state, 

where A" = ?ik, ''^k = cj^Ck is the total particle density, the average is taken in 
|BCS,^), and ^k = (^k) = I'^^kP- Here we have used — rik, implying (rik) = (j^k), 
as well as {hqUk) — {nq){hk) for k 7^ q due to the product nature of the BCS state. 
Since < ^k < 1 and rik = only for a non-extensive set of modes, the scaling of the 
numerator is ~ A^, and that of the denominator ~ A^^ in the thermodynamic limit. It 
is the mean-field (product) nature of the state which is responsible for this scaling. 

The fixed phase wave function is an exact solution of the number conserving 
equation (if the particle number is not fixed to a specific A^), since |BCS,6') oc 
exp(Q;G^)|vac) = ^^(Q;G''")^/A^!|vac). The above discussion shows that for fixed 
A" — > 00, the fixed phase wave function approximates the fixed number one in a 
controlled way. We will rely on this fact, which is solely related to large A^, in the 
mean- field theory to be discussed next. 
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Appendix A. 2. Mean-field theory 

Here we describe a mean-field theory for quartic, number conserving fermionic master 
equations. The mean-field description, represented by an effective quadratic master 
equation, is valid in the thermodynamic plus long-time limit, where the system is already 
close to the steady state. It is instructive to compare this mean-field theory to BCS 
theory for Hamiltonian ground states. The role of the ground state is played by the 
dark state. The long time limit is analogous to a low energy limit in condensed matter 
systems, where the description of interacting fermions in terms of quadratic effective 
BCS Hamiltonians becomes appropriate. We emphasize that unlike the BCS problem of 
locally interacting fermions, our mean-field theory is controlled by the fact that the exact 
dark state for the interacting quartic dissipative problem with fixed particle number is 
known. The effective linear Lindblad operators take the form of annihilation operators 
for the dark state, in complete analogy to the Bogoliubov operators for the BCS problem. 
The momentum dependent rate premultiplying the Lindblad form describes the damping 
in the vicinity of the dark state, replacing the excitation eigenenergies in the BCS 
Hamiltonian. The mean-field theory allows us to calculate the damping rates and 
number equation from the microscopic quartic description. 

In our mean-field approach, as in standard BCS theory we give up exact particle 
number conservation and work in the fixed phase ensemble, exploiting the discussion 
above. Our mean-field ansatz is then defined by a factorization of the density matrix 
in momentum space, p — HkPk, where pk describes the mode pair {k, — k}. This is 
motivated by the form of the coherent state steady density matrix po — |BCS, 9) (BCS, 9\ 
sharing this property. We implement this ansatz, focusing e.g. on the recycling term, 

^[p] 3 ^^kPA = XI ^qi'"q2<3'^q4 4i"q2^43«q4<^(qi + Qs - q2 - q4) (A.13) 
k qi,-,q4 

= X VpUq^u*^^vlalaq^pal^ap5(qs - 

q2,q37^p 

+ VqiUpU*pV*^^al^^appalaqJ{qi - q^) 

+ VpUq^u*_pV*^^ala^^palpaqJ{ci2 + q.4} 

q2,q47^p 

+ VqiUpU*^y_palfypal^^a_p6{cii + qg) + {p ^ -p} + h. o. t. 

qi-qsT^p 

For the last equality, we have selected the contributions which are of second order in the 
operators for some selected momentum mode p and neglect the higher order terms in 
the following, as appropriate in the thermodynamic limit. Now we perform the partial 
traces using [pk,Pq] = [Pk, Oq] = [Pk^o^ = for q 7^ k, since pk contains an even 
number of fermions. We furthermore reorder the fermion operators using that the sums 
do not contain the modes ±p; i.e., we only have to use the anticommutation property 
{cq, = for q 7^ p. That is, the procedure will only yield signs depending on how 
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many reorderings are necessary: 

Tr^p>C[p] 3 Tr^plVpVlalppttp WqWqaq JJ + UpU^ttpPpal ^ f qf qfl^ JJ PkOq 

Qt^P kT^p q^p k^p 

- VpU*_palppalp ^ liqi;* qCq JJ PkO-q 

q^p k^p 

- UpV*_papPpa^p v_qaLq JJ pkO^ + {p ^ -p}] 

q^p k^p 
= t^p^paJ,Ppap ^q^q(«i«q) + '^p'^p^pPp^^l Yl ^^K^^^^l^ 

<i¥=p <i¥=p 
- VpU*_palppalp Y Wq^^-q(a-qaq) - UpV*_papPpa_p Y «q^-q(«ia-q) 

+ {p^-p}, (A.14) 

where in the last equahty we have chosen the normahzation TrkPk = 1 and the cychc 
invariance of the trace. In the thermodynamic hmit, X]q-^p(-) = I^q(-)- Close to 
the steady state we can evaluate the correlation functions using the stationary form 
equation (A. 10) for the state, which take a factorized form, i.e., (aqOq) = l^qP; (^q^q) — 
l^qp, (a-qOq) = —u'^Vq, (OqoLq) = —v^Uq. Wc cau thus factorizc a real positive 
number kq from the last equation, such that mean-field Liouvillian operator for the 
mode pair {k, — k} takes the form 

CkW] = Tr^p-^^b] ^ i^oY (^<^kPk^ik - ^{^Ik^^k, Pk}) (A.15) 

a=± ^ 
= Yj '«k(-t'(TkPk-t'l.k - ;^{-^ak-^(^k,Pk}), 



q q 'I 

with -Lk,Zk specified in Eqs. (A. 5), (A. 8). 

In addition to the effective dissipative rate, we can also calculate the amplitude of 
the relative coefficient a in dependence of the particle filling n — X^q(aqOtq) from mean- 
field theory. While for the number conserving equation the particle number operator 
is a constant of motion due to [^j,iV] = (i.e., the expectation values of any power of 
the total particle number arc conserved, dt{N"'') = 0), no such property is found for 
the mean-field dynamics. However, for any initial average particle filling n = {N)/L'^ 
there is an a = re'^ such that the first moment, i.e., the average particle density, is 
conserved, and dt{N) — holds. In other words, the value of a fixes the average particle 
number, which is intuitive since a measures the relative strength of creating fermions 
vs. annihilating fermions. The number equation of state for which n{t) — const, is given 

by 

- = E l^.l' = E rJH^- (A.16) 



q q 



|Vq|"^ + \aVq\ 
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While, therefore, the modulus |q;| = r is fixed by an additional physical condition from 
the initial state, no such condition exists which fixes the phase 9. This is physically 
meaningful, since the original number conserving operators are invariant under global 
U{1) phase rotations. 

The interpretation of a = re^^ therefore is the following: (i) the macroscopic phase is 
fixed spontaneously by the dynamics, i.e., it indicates spontaneous symmetry breaking of 
U (1) phase rotations in the dissipative setting, (ii) The amplitude relates to the average 
particle number in the sample, which can be chosen such as to match the average particle 
number of the number conserving setting n — N/ in the initial state. 

Appendix A. 3. Explicit derivation of the general form of Majorana zero-damping 
modes for translation-invariant dissipative processes 

In this section, we provide an explicit derivation of equation (54) in the framework 
introduced in section 5.2 (using the same notations). 

An arbitrary Majorana operator 7 = qqj^ ]-,g expressed, in terms of the fermionic 
creation and annihilation operators and Oj corresponding to distinct lattice sites i of 
the system, in the form 

7 = ^^(ctjOj + /i.e.), (A-1'^) 

j 

where the sum runs over all sites j belonging to the system and aj G C with |«j p = 1 
(so that {7,7^} = 2 as required for a Majorana operator). The operator 7 then 
corresponds to a Majorana zero-damping mode of the dissipative dynamics described by 
the translation-invariant Lindblad operators Lj = ^^gj^j) Wj-i aj+Vj^i Oj (equation (52)) 
if and only if the following conditions are satisfied (see section 2.7): 

= {Li,7} = {L„7,} (foralH), (A.18) 

where 7^ denotes the restriction of the Majorana operator 7 to the sites j G X(i) where Lj 
acts non-trivially. These anticommutation conditions obviously constrain the possible 
form of 7i in a way that solely depends on the form of Lj. Since the form Lindblad 
operator Lj does not depend on i due to translation invariance, 7^ must be of the form 

7i = XI iPj-i'^iO'j + ^-C-)' (A. 19) 

where (3j-i are complex factors {(3j-i being a shorthand notation for /3(rj — r,), as defined 
in section 5.2) which, for consistency with equation (A. 17), must satisfy the following 
set of equations: 

/3j_iai = aj] = (A.20) 

(3i_jaj = = "I, (A.21) 

leading to a single "consistency relation" 

/3k-i = /3k-jPj-i, (A.22) 
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which is vahd for any triple of indices k) associated with lattice sites belonging to 
the system, with = /3o — I. Since /3j^i — {/3i^j)~^ and = the factors 

must be real. 

Equation (A. 22) implies that any Majorana zero-damping mode 7 can be fully 
constructed from the only knowledge of d real factors /3e„, where e„ (n = 1, 2, . . . , 0?) are 
primitive vectors associated with the d-dimensional Bravais lattice on which the system 
is defined. Such a construction can be made starting from any lattice site i, choosing a 
phase (f)i e [0, 27r) for the coefficient ai of 7 (see equation (A. 17)) corresponding to that 
particular site. This allows us to express 7 in the generic form presented in the main 
text (see section 5.2): 

d 

J = J2 {f3e,r{Pe,r ■ ■ ■ {PeaTM^i + m„e„) + h.c. ,(A.23) 

{mx,m2, n=l 

where A/" > is a normalization factor and {mi, 7712, . . . ,^^1 a set of integers defined 
such that the vectors + Ylt.=i "^n^n span the positions of all sites in the system *. 

The explicit form of the elementary "building blocks" /3e„ defining the form of 7 
can be found by expressing equation (A. 18) exphcitly using equation (A. 19), namely, 

0= ^ ^ ^Uj^iaj^Vj_ia],l3k-iOiiai + I3k-ia*a}i^ 

kex(i) 

^ ] {(^i '^j—iPj—i ~l~ (^iVj—iPj—i) 

iex{i) 

where in the last step wc have symmetrized the expression using the fact that X{i) 
contains pairs of sites located symmetrically around the center of symmetry Si 

of Li (see section 5.2). The fact that the Lindblad operators are symmetric around Si 
and that the dissipative dynamics leads to a pure steady state as assumed in the main 
text implies that Uj-i — —u^-i and Vj-i — Vj_i (see also equation (39) and discussion 
thereof). Equation (A. 24) therefore reduces to 

= J2 - Pj-i) + vj-i{Pj-i + /3j_,)] , (A.25) 

which is the result presented in the main text (see equation (54)). We note that this 
equation can be further simplified when the center of symmetry Si of Li coincides with 
a lattice site. In that case, = = (/3j_j)~^ and we find 

0-Vo + IJ2 [^'"^'^J-^{PJ-^ - iP^-i)-') + + (A.26) 

iex(i) 

We remark that equation (A.25) (as well as equation (A.26)) can be expressed in terms 
of the factors only using equation (A. 22). If there exists a set {/3e„}n=i,2,...,d of 

*Note that 7 = 7^3 defined up to a sign, which allows us to introduce the phase 0j as e''^'/^. 
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real factors ^e„ satisfying equation (A. 25) for some particular phase 0, G [0,27r), 
the system can support at least one Majorana zero-damping mode of the form (53). 
Otherwise there can be no such mode. In order to assess the possibility of having 
Majorana zero-damping modes, one therefore has to solve a polynomial equation in d 
real variables {{(3e„}n=i,2,...,d) with coefficients that may be real or complex depending 
on the coefficients e~''^^Uj-i, since vj-i can always be chosen as real (see section 4). 
The existence of solutions thus depends on the phase 0j that is chosen as an "initial 
condition" for the construction of 7 starting from site i according to equation (A. 23). 

We finally discuss the physical meaning of the phase (pi in the two types of dissipative 
systems that are most relevant to this work (see section 4): (i) ID dissipative systems 
with TRS (symmetry class BDl) and (ii) 2D dissipative systems without TRS (symmetry 
class D) *. In the ID case, TRS constrains the coefficients Uj_i to be real up to a global 
phase and the phase (pi crucially allows to compensate for such a phase so that all 
coefficients e"'*^* real. Equation (A. 25) then reduces to a univariate polynomial 

equation with real coefficients which, as opposed to a univariate polynomial equation 
with complex coefficients (as would generically be obtained in the absence of TRS), 
can have robust solutions (i.e., solutions that survive if the coefficients Uj^i and Vj^i 
are slightly modified). In the 2D case, Vj-i typically exhibits an isotropic (or s-wave 
symmetric) form due to physical constraints (see section 4) and TRS can only be 
broken if Uj^i = u{rj — r^) transforms under rotations as an eigenfunction of angular 
momentum with odd integer eigenvalue £ ^ 0, namely, ^(r) = u{r,Lp) ~ /(r)c^'^''' with 
polar coordinates {r,ip), f{r) being an arbitrary positive function of r. We then find 
e~^'^*Uj-i — e~^'^^u{r,ip) = u{r,ip + (pi/i), showing that modifying the phase (pi simply 
corresponds to rotating the reference frame or, equivalently, to modifying the orientation 
of the edge of the system. If there exists a solution of equation (54) for some 0j, there 
must exist a solution for any (pi e [0, 27r) or, equivalently, for any orientation of the 
edge. TRS breaking thus crucially allows for the existence of Majorana zero-damping 
modes at edges that are curved arbitrarily, as expected physically. 
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